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ABSTRACT 


A theoretical engineering thermal analysis of man-made "permafrost 
islands" as offshore drill sites is conducted by developing a two- 
dimensional numerical model from the heat conduction equation. The 
model is formulated by solving an equivalent Pane principle 
based on the finite element method. Phase transformation due to 
‘pene ios and thawing is assumed to occur isothermally to promote 
Solution economy. Solutions in the time domain are obtained by 
the Crank-Nicolson method. 

Boundary conditions may be arbitrarily imposed or simulated at the 
air-surface interface. Surface boundary conditions are simulated by 
ate a one-dimensional surface energy exchange model which is coupled 
to the two-dimensional model. Local meteorological data is employed 
for the surface energy simulation. Successful tests were conducted 
that compared the numerical model solutions with exact solutions to 
freezing front propagation. Also, comparisons were made between model 
' results and in situ measurements of soil temperatures at Fairbanks, 
Alaska. Finally, the model was applied to the problem of offshore 


"permafrost island" drill sites. 
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Chapter 1 


INTRODUCTION 


1.1 Alternatives to Conventional Offshore Drilling Structures 


Recent studies indicate that offshore regions in the Beaufort and 
Chukchi Seas may contain extensive oil and gas reserves. Some geologists 
Beiteve that offshore areas tend to be more productive than adjacent 
land areas. If this is true for arctic waters, the implications are 
Staggering considering the recent petroleum discoveries along the arctic 
coast. It can be expected that exploratory offshore drilling in the 
Beaufort Sea will commence sometime within the forseeable future with 
initial probes restricted close to shore. 

OF fshore drilling experience in hostile arctic environments has 
been limited to Cook Inlet, Alaska. Although annual growth of sea ice 
in Cook Inlet is only 3 feet, design wind, wave, and earthquake forces 
on offshore structures have been found to be negligible when compared 
with expected lateral ice forces. Obviously, the coastal environment 
in the Beaufort Sea is much more hazardous to ocean structures since 
annual sea ice reaches thicknesses of 7 to 9 feet and is subject to 
continuous motion due to the driving forces of wind and current. 

Other dangers include the possibility of pressure ridges, ice islands, 
or multi year ice striking a structure. 

Assuming an adequate knowledge of critical design loading due to 
ice, the application of known engineering techniques for offshore drill- 
ing structures could probably be physically implemented but, at an 


economically prohibitive expense. The arctic thermal regime may permit 
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consideration of two additional alternatives for development of offshore 


drill sites: 

1. Construction of man-made permafrost islands from dredged sea- 
floor materials that are frozen into an icy matrix strong 
enough to resist lateral forces exerted by annual sea ice 
movement. 


2. Construction of structurally reinforced man-made ice islands 


grounded at the desired location. 


1.2 Scope and Objective 


This report is addressed to the first alternative. The primary 
problem, thermal analysis of nearshore man-made permafrost islands, is 
considered. The objective is to develop a versatile engineering tool 
_ based on known thermodynamic theory and make meroeet ntd assumptions to 
facilitate engineering application. A mathematical model is formulated 
for computer solution utilizing the finite element method. Numerical 
solutions are compared with analytical solutions prior to application 
to example permafrost island problems. This report is restricted to a 
theoretical study relying on available data for soil thermal properties 
and available historical records for meteorological conditions. The 
results of several long-term simulations are presented. They are 
considered to be of value to engineers desiring design information on 
equilibrium freezing front propagation in various soil structures 


subjected to offshore arctic meteorological conditions. 
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Chapter 2 


FINITE ELEMENT FORMULATION OF TRANSIENT HEAT CONDUCTION WITH PHASE CHANGE 


2.1 Freezing and Thawing of Soils 

Heat transfer with phase change due to freezing or thawing is the 
principal thermal process of interest in applied arctic soil engineering 
problems. When attempting to mathematically quantify or simulate the 
physics of these processes, the engineer must assess what level of model 
sophistication is necessary in the context of the proposed application. 
The primary factors affecting his decision are the extent to which the 
physical properties of the soil system are known and the quality of 
information available for determining external boundary conditions. 

As a general rule, the physical parameters of the soil system will 
dictate how the phase change problem should be solved. Three alterna- 
tive methods are considered for modeling the latent heat of fusion 
problem for saturated soils: 

A. Phenomenonological 

‘B. Thermodynamic 

C. Analytical - isothermal 


Each method and its inherent limitations is briefly discussed. 


2.1.1 Phenomenonological Modeling 

It is well known (Keune and Hoekstra, 1967) that certain moist or 
saturated soils do not undergo isothermal phase change. Lukianov and 
Golovko (1957) proposed a phenomenonological equation to describe such 


a process in a soil-water-ice system: 


a 
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a9 ao 
Ox Kx) = For (1) 


where x = distance from the air-soil interface (ft) 
t = time (day) 
> = temperature (°F) 
K = K (x,t,o) = thermal conductivity (BTU/ft-day-°F) 
F = apparent volumetric heat capacity (BTU/£t>-°F) = C - L dG/d¢ 


Cy = C, (x,t,6) = volumetric heat capacity of the soil-water-ice 
system 


G(?) = Ice content (1b/£t>) 


L = latent heat of fusion (BTU/1b) 


Equation (1) mathematically describes the freezing-thawing phenome- 
non in a one-dimensional zone of finite thickness. The freezing of a 
_small differential element is accompanied by a temperature dependent 
release of latent heat which decreases the time rate of change of 
internal energy if net heat extraction from the element remains invariant. 
The phenomenonological approach is equivalent to introducing an "apparent 
volumetric heat capacity" as defined above. 

Although phenomenonological modeling of phase change remains 
analytically intractable, it is possible to numerically solve equation (1). 
Nakano and Brown (1972) developed a one-dimensional finite difference 
scheme from equation (1) and satisfactorily simulated heat conduction 
in tundra soils near Barrow, Alaska. It is important to stress that 
the success of any phenomenonological model is heavily dependent on 


the following factors: 





1. Unfrozen soil water content as a function of temperature must 
be known or experimentally determined. 

2. Numerical simulation of phase change requires small time incre- 
ments to avoid solution instabilities caused by the injection 


or removal of latent heats. 


3. Applied boundary conditions must be known on a micro-scale. 


2.1.2 Thermodynamic Modeling 

Thermodynamic modeling of soil phase change also incorporates 
equation (1) to predict soil column temperatures as a function of time 
but does not rely on experimentally determined soil freezing curves to 
simulate the temperature dependent evolution of latent heat. Recent 
research (Koopmans and Miller, 1966) has shown that many soil moisture 
‘characteristic curves and soil freezing curves are similar. Since a 
soil moisture characteristic curve defines the relationship between 
water content and soil pore water pressure, detailed calorimetric 
analyses to determine unfrozen water content in soils as a function of 
temperature can be avoided if temperature-pore water pressure relation- 
ships during phase change are known. 

It is possible to calculate soil water pore pressure in soils under- 
going phase transition by assuming that pure ice and water can only 
co-exist in thermodynamic equilibrium as defined by the Clausius - 
Clapeyron equation for a bulk fluid (Sears, 1964) 


LdT 


ap= TQiav) @) 
where L = latent heat of fusion 
T = existing absolute temperature 
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ao 


ee 


v'' = specific volume of ice 


specific volume of water 


< 
u 


incremental change in pressure 


Q 
a] 
u 


It is sometimes assumed that the Clausius - Clapeyron equation applies 
to unsaturated soils and that ice pressures in the soil pores are at 
atmospheric pressure. It is probable that this equation does not apply 
to saturated soils without modification. 

A typical numerical procedure for thermodynamic simulation of soil 
freezing is to compute a soil temperature profile from equation (1), use 
this to obtain soil water pore pressure by numerically integrating 
equation (2), and enter the soil moisture characteristics curve to 


determine the change in water content. The net change in water content 


. over a given time increment determines the magnitude of latent heat 


that must be released in the following time step. A new or updated 
apparent volumetric heat capacity is computed to continue the temporal 
solution of equation (1). This method of discretely modeling phase 
transformation requires the utilization of small time increments to 
avoid unstable oscillations in computed temperatures. Furthermore, 


soil moisture characteristic curves must be experimentally determined. 


2.1.3 Analytical Modeling of Isothermal Phase Change 


Analytical solutions commonly employed in the engineering analysis 
of phase change problems are simplified versions of the theory of ice 
growth introduced by Franz Neumann in the 1860's. Examples are Stefan's 
(1889) and Berggren's (1943) applications of the Neumann theory to ice 


growth and frost penetration problems. Since any theoretical formulation 
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of phase change (commonly referred to as the "problem of Stefan") is 
nonlinear, classical solution techniques are not applicable (Carslaw 
and Jaegar, 1959). 

Because of the mathematical complexities involved, the practicing 
engineer is usually forced to seek a published analytical (exact) 
solution that represents the best approximation to his physical problen. 
A significant disadvantage of mathematically exact solutions is the 
limitation imposed by simplifying assumptions and solution space geo- 
metry. To analytically solve the isothermal phase change problem, it 
is usually necessary to assume a semi-infinite isotropic medium initially 
at a uniform temperature (i.e., simplified boundary conditions are 
chosen so that the governing equation of state can be solved). The 
imposed boundary condition is a sudden change in temperature that subse- 
quently induces phase change in the medium. It is the restrictive 
imposition of simple boundary conditions that severely limits the 
general application of exact solutions. 

Several Russian investigators have been successful in analytically 
modeling more realistic boundary conditions. Kolesnikov (1946) presented 
a derivation which incorporated the major hydrologic energy exchanges 
affecting the growth of sea ice, and more recently, Portnov (1962) 
developed an exact solution for isothermal phase change with time 
dependent boundary temperatures. 

Because of the relative paucity of exact solutions and the complexi- 
ty of modern engineering problems, most research efforts have been 
directed toward achieving approximate methods of solution. The most 


effective have been finite difference and finite element methods which 

















<a 4 3 7 j 
= = af ¢ 
a. Bcc . ihe, ae”) 
‘glo? Fp wolddrq” ois ne ot be re ib 
nal id ; 41. ws o 


a 


ws . qqs fon ots sbupiadtesy nity Ree 






} \Wbevloyal-<citixotqaco lazltnsedsem ang Mod 







j dearest Frerfe bid\ ae « Neo? o3 heart Vi Dewars be 


al 
TF 25 


7 


q aia Of Molicainerqe? 2250 of} tenet 


a fa. Tonka “fiaolopmadszie Yo one tevbneake " 
Meksibor brs ancizqevees. gilyitignte faa 
(ag nda Jaerrndsoet uds svFou ssclieall 


shegheg olqoxiozi eiinital-imse os anijgzs oF (iuanesen 
ey eis =o LT bine» ytehneed Loitilqeie , 9.2) sceegeed aa 


++, ree “a 
4 Lyr oat iclict kas O7B7E io AG (rsovog ont taal! 
; : ~~ 

tue seu walnkseqnes mi eginato libbuz ‘ 


i ‘evizoFaaast B02 31) mE Bee “Sir ot vxnans shane = 






mi 














us 
‘ 
*” 


, 
7 


“iti fred sik 


¢ 
& ont ddmht “Sexover S5:ft tnottiiacd atin elQake 
, loz iorte’ 3 wiveadd 
‘Nbviottinas Ricfwiz2eoone med over etotsgisesvini agivaull . ot ave 
Bssnekorg {ShOl) vortices !os andi> thiios . ¢rofnucd cl serine sin 
tepusiiare yore uivoiorivi sofas 49; e305 iu falew trod * ri 
7 


(Seet} Plats ie Cail’ aves [> se 


edit 2 Nk eucods erbigq Laurmiistocr : c5is iE: 
‘ 
a Poe Atay I. bo® Ag 
7 A , 
173e@ SVAN @f To » 4 | it] 
a Oc ( i 
a Jit fan é it¢ \ ‘ f 4 





— 


rely on analytical solutions as a means of testing numerical techniques 
prior to simulating physical problems that remain analytically insoluble. 
Hence the primary importance of known theoretical solutions is that of 
investigating the accuracy of a numerical procedure whose validity has 


not been theoretically established. 


2.2 Derivation of a Two-Dimensional Finite Element Model 

An approximate numerical solution to the parabolic differential 
equation of transient heat conduction in two dimensions is developed by 
employing the finite element technique based on the so-called Rayleigh- 
Ritz formulation. Phase change is assumed to occur isothermally to 
avoid previously discussed difficulties common to thermodynamic and 


phenomenonological methods of solving the latent heat problem. 


2.2.1 Variational Finite Element Formulation 

The variational functional approach for the finite element formula- 
tion of transient heat conduction relies principally upon the theoretical 
contributions of Biot (1970) who introduced a unified variational 
analysis of heat transfer analogous to the energy theorems of classical 
mechanics. Given the partial differential equation of two-dimensional 


transient heat conduction in Cartesian coordinates 
) od e) dd ee 
3x Kay? yi ty Yy3y? + G pcr (3) 
where 4 = temperature 


Ky and are the thermal conductivities in the x and y 
direttions, respectively 


G = G(x,y,t), a variable heat source or sink 


pe = volumetric heat capacity 
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t = time 
the transition to an equivalent variational statement is sought by 
applying an extension of the fundamental lemma of variational calculus 
(Weinstock, 1952). The theorem is applied by noting that if in a given 


time interval, the functional 


= [SEQ 9g, 36: 
x Rf GY bose addy (4) 


is to be minimized over a bounded region R, a necessary and sufficient 
condition is that at any given instant of time the unknown function 
$(x,y,t) satisfy the Euler equation 


Eero) + x oem) - oe = 0 (5) 
Ox ‘0 (06/7 0x): dy 0 (00/ay) oo 


subject to $ obeying the same boundary conditions of equation (3). The 
problem is to obtain a functional which when substituted into equation 
(5) reduces to equation (3). 

Assuming se remains invariant, an equivalent variational formula- 
tion of equation (3) is obtained if the functional 


x = FLf (Go)? + Ky (Go)? + 2(oege - G)¢}dxdy (6) 


can be minimized for all admissible states of >. For equation (6) to 

be a true minimizing principle, the second variation of the functional 
with respect to the unknown state variable must be greater than zero. It 
can be shown that this is indeed true for transient heat conduction 
(Myers, 1971). It is less restrictive to require the functional to be 
stationary (i.e., only the first variation is zero) and demonstrate 


that the numerical solution will converge to an exact solution. 
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Equation (6) is formally solved by discretizing the domain R into 
numerous subregions or "elements" and assuming that the temperature 


th 


distribution within each element can be approximated by an n~ order 


polynomial. In this report, a linear polynomial is selected to 
approximate the state variable in a two-dimensional finite element mesh 
composed of arbitrarily shaped triangular elements. 

Upon substitution of a linear temperature distribution shape 
function into equation (6), the first variation of the functional is 
evaluated by applying Leibnitz's rule for differentiation of integrals. 
Thermal properties are assumed to be constant within each element although 
abrupt changes in properties between elements are permitted. Summation 


of the integrated contributions of all elements produces a system of 


simultaneous linear first order differential equations of the form 


[k]{o} + [c] (2% = {c} (7) 
where [K] = system conduction matrix 
(C] = system capacitance matrix 
{G} = system load column matrix 
{o} = nodal temperature column matrix 


The system conduction and capacitance matrices are banded and symmetric, 


Significantly reducing computer memory requirements for storing arrays. 


2.2.2 Incorporation of Isothermal Phase Change 


Simulation of heat conduction and phase transformation in two 
dimensions involves a minor modification of the variational formulation 
of transient heat conduction. As previously discussed, the solution 


domain is discretized by a finite element mesh composed of triangular 
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_ elements. The mesh geometry is stored in computer memory by reading in 
nodal points, their respective global coordinates, thermal properties 
of each element (frozen and unfrozen), and moisture content of each 
element. After this has been accomplished, numerical simulation of 
phase change may commence. os 

Consider a common node point "A" in a segment of a two-dimensional 
finite element mesh (Figure 1) whose temperature is dependent upon any 
arbitrary boundary condition imposed along S. Phase change is not 
permitted to occur at node "A'' until a requisite amount of latent heat 
has been exhausted or added. The prevailing quantity of latent heat 
for node "'A'' is determined by adding one third the weight of material 
(subject to phase transformation) of each contributing element to a 
nodal latent heat accumulator array assigned to node "A". 

As an example, consider the cooling and freezing of node "A". 
After each time increment that a solution is being generated, the temp- 
erature of node '"'A'' is scanned to determine whether or not it has 
dropped below a prescribed point where phase transformation is assumed 
to occur (TFAZ). If the node has not been previously frozen and the 
computed temperature has dropped below TFAZ, then the quantity of 
latent heat evolved during the previously computed time step becomes: 

DELQ = C,, (TFAZ - Computed Temperature) * Volume (8) 
where Gy = a weighted volumetric heat capacity based on the heat 

capacities of the soil-water mixture assigned to node "A" 
that must undergo phase change. 
Volume = volume of soil-water assigned to node "A" 
The quantity DELQ is subtracted from the latent heat accumulator for 


node "A'' and if more latent heat must still be exhausted, the computed 
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temperature at the node is brought back to TFAZ prior to obtaining a 
new solution at the end of the next time increment. This is analogous 
to simulating phase change as an isothermal process. The time rate of 
latent heat generation is governed solely by the solution of the heat 
conduction equation and remains insensitive to large time steps. 

If the required amount or more than the required amount of latent 
heat has been extracted, node "A" and the volume of soil-water assigned 
to it are considered to be frozen. The thermal properties of the 
elements common to node "A" are then updated based on the volumetric 
proportions of soil-water and soil-ice present in the respective 
elements. If an excess amount of latent heat is removed, the residual 
is used to calculate how far below TFAZ the temperature at node "A" 
should be. A weighted C¢ is used for this computation. All weighted 
volumetric heat capacities (Cy and Ce) are computed from the thermal 
properties of elements common to nodal points undergoing phase change. 

The thawing process is handled in a similar manner. Provisions are 
made to monitor which nodes are frozen or thawed so as to avoid re- 
freezing frozen nodes. Boundary nodes remain unaffected by the checking 
routines that scan nodal temperatures. After phase transformation at 
a node occurs, the latent heat accumulator assigned to that node is 
recomputed. Recomputation permits the simulation of long-term cyclic 


freezing and thawing. 


2.2.3 Imposition of Boundary Conditions 
The variational functional expressed by equation (6) automatically 


incorporates boundary conditions of either specified temperature or 
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zero heat flux. Inclusion of additional modes of heat transfer that 
may occur over various portions of the boundary surface, shown in 
Figure (1), is accomplished by adding appropriate surface integrals 

to equation (6) which upon minimization yield the desired boundary 
conditions. The surface integrals for specified heat flux, convection, 


and radiation heat transfer across the boundary are 


f adds + Lf K(4? - 2bn6)ds + of (Lego® - egodoyds (9) 
Be 2B. Be > 


where q = specified heat flux across the boundary 
@ = boundary node temperature 
Doo = sablonbiceipetnenris near the boundary 
K = convective heat transfer coefficient 
o = Stefan-Boltzmann constant 
€p = Emissivity of the boundary surface 
E, = Emissivity of the ambient surroundings 
Outgoing heat fluxes are considered positive. 

The combination of equations (6) and (9) and their subsequent 
minimization produce a generalized equivalent variational statement 
describing transient heat conduction with its associated boundary 
conditions. Substantiating mathematical details are provided in 
Appendix A. Minimization of surface boundary integrals produces terms 
that are combined with the system conduction and load matrices of 
equation (7). Nonlinear radiation heat transfer terms give rise to 
solution difficulties since unknown temperatures are introduced into the 


System matrices. Methods of circumventing this problem include lineari- 
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zation or modeling net radiation exchange as a specified heat flux. 


These methods are discussed in greater detail in section 3.2.2. 


2.2.4 Solution in the Time Domain 
Linearization or removal of nonlinear terms from the system 
matrices allows preservation of the previously discussed system of 


ordinary differential equations represented by 
[k]{$} + [c] (2% = {ch (20) 
ot 


and greatly simplifies the time domain solution. After a considerable 
amount of experimentation with various numerical time domain solution 
techniques, the Crank-Nicolson method was selected since it provided 
the most satisfactory results in terms of optimizing time step sizes, 
of providing acceptable solution accuracy, and of minimizing numerically 
induced oscillations. 

The Crank-Nicolson method utilizes the arithmetic mean of the 
derivative of the state variable at the beginning and end of each time 


step to move the solution ahead in time. Thus 


{gyV*} = 9)” + SEctonY + (6F"*4) (11) 


where At = time step size 


v = a reference point in time 
Upon premultiplication by the capacitance matrix and substitution of 


equation (10), equation (11) becomes 
CIK] + @1c]) {6}? = GZIc] - [Kk] fo” + 2{ch (12) 


Since {o}” is the initial or previously computed temperature 


distribution and the system matrices [K], [C], and {G} remain constant 
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in each time step, equation (12) can be further simplified to the 
matrix equation 

(w{o}’*! = {2} (13) 
where [W] = a symmetric matrix in band form. 
Instead of employing a time consuming matrix inversion algorithm, the 
solutions, {¢}¥*!, are computed by the Gaussian elimination method which 
exploits the symmetric and banded structure of the [W] matrix. 
2.3 Comparison of Finite Element Solutions with Analytical Solutions 

Involving Phase Change 

Because the method formulated to solve the latent heat problem is 
a numerical approximation, comparison of finite element solutions with 
several exact solutions to simple isothermal phase change problems is 
justifiable. Two problems are considered: one in which thermal para- 
“meters remain unaffected by phase change and secondly, one where phase 
transformation produces a significant variation between frozen and 
unfrozen thermal properties. 

.2.3.1 Simulation of Isothermal Phase Change with No Change in 

Thermal Parameters 

As a first example, the problem of the solidification of a homo- 
geneous slab of finite thickness is considered. An exact solution 
discussed by Allen and Severn (1952) is used for comparison purposes. 
In Figure 2, a two-dimensional finite element mesh for analyzing the 
problem is shown. Figures 3 and 4 summarize the finite element results 
obtained from two different time domain solutions: implicit and Crank- 
Nicolson. Good agreement between both numerical solutions and the exact 


solution is evident. Note that the Crank-Nicolson time domain solution, 
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Figure 2. 


Figures 3 and 4. 
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Solidification of a slab of liquid: finite element 
mesh. 


Distance between nodes 1. and 19 = L = 60. AX = L/6 = 
K = 1; pe = Y2 (thermal properties remain constant) 
Uniform initial temperature = 45° 

Temperature of phase transformation = 45° 

Latent heat of transformation = 70.26 pc/unit volume 


The entire perimeter with the exception X = L is 
thermally insulated. 


The end X = L at t = 0* is dropped to and subse- 
quently maintained at a temperature of 0°. 
Solidification of a slab of liquid: progress of 


freezing front through the slab. 


Inset diagram illustrates nomenclature in each 
figure. 


The freezing front is denoted by line F at a typical 
time t, and the frozen zone is shown shaded. 


Exact solution for the equation of the solidification 
locus is: 


ee 
(X-L)? = SF(1.07)t 


as obtained from Jones' theory (Allen and Severn, 
1952) 
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Figure 4, permits similar solution accuracy with time steps twice as 

ae as those employed in the implicit time domain solution, Figure 3. 

This example demonstrates the accuracy of finite element discretiza- 
tion of the freezing process. Judicious selection of time steps and 
finer mesh subdivision will improve solution accuracy with increased 
computer simulation time. A much more rigorous and realistic test is 
to consider the following theoretical problem where substantial changes 
in thermal properties of the conducting medium occur upon phase change. 

2.3.2 Simulation of the Neumann Theory of Freezing Front Propaga- 

tion in a Saturated Soil Column 

The application of the Neumann theory to the prediction of frost 
penetration depth in soils is well known. A transcendental equation is 
solved to obtain a constant of proportionality whose magnitude is 
dependent upon the thermal properties (frozen and unfrozen) of a semi- 
infinite soil column, initial uniform temperature distribution, and 
imposed initial constant boundary temperature. This constant yields the 
exact solution for freezing front propagation when multiplied by the 
square root of elapsed time. 

Comparison of a one-dimensional finite element simulation with a 
Neumann solution for freezing front propagation is shown in Figure 5S. 

_ The Neumann solution was extracted from an example cited by Jumikis 
(1966). The freezing front locus is plotted versus the square root of 
time to demonstrate the stability of the finite element solution and 
enable computation of the numerical constant of proportionality. The 
theoretical constant of proportionality is .575 £t/day~1/2 and the finite 


element solution constant of proportionality is .56l ft-dayl/2, With 
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Comparison of a variational and analytical solution 
to the Neumann problem. 


A. 


B. 


Physical constants: saturated sand at a void 
ratio of .50; specific gravity = 2.65. 


Thermal parameters: 
1. Unfrozen: 

K, = 25.7 BTU/ft-°F-day 

C, = 42.7 BTU/ft?°F 
2. Frozen: 

Ke = 32.2 BTU/ft-°F-day 

Ce = 29.3 BTU/ft=°F 
Initial and imposed boundary conditions: 
Initial uniform temperature distribution of 
36°F. A constant boundary temperature of 14°F 
is maintained. Isothermal phase transformation 
is assumed to occur at 32°F. 
Analytical and variational solutions are: 
x = .575 vt (theoretical) 


x = .561 Yt (finite element) 


where x = frost penetration depth (feet) 


t = time (days) 

Variational solution was obtained by employing 
a one-dimensional finite element mesh. Spacing 
between nodes was one foot with numerical 
computations conducted in .5 day increments. 
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such a small discrepancy in constants of proportionality, it would 

take an elapsed time of 15 years before the error between computed and 
theoretical frost penetration depth exceeded one foot. At that time the 
theoretical freezing front would have penetrated 42.6 feet below the 


ground surface. 
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Chapter 3 


SIMULATING AIR-GROUND INTERFACE ENERGY EXCHANGES 


3.1 Methods of Approximating Heat Transfer at the Air-Ground Interface 

The thermal processes occurring at the interface between a soil 

system and the atmosphere are often too complex for an exact descrip- 
tion; yet, their very existence is crucial to conducting an analysis 

of the thermal regime. The dominating processes are absorbed incoming 
solar radiation, net long wave radiation exchange between the atmosphere 
and surface boundary, turbulent heat transfer, evaporation, snow and 
vegetative cover effects, and the soil temperature gradient at the air- 
ground interface. There are three alternative methods that may be 
employed to simulate the net effect of these processes: 

1. Impose a soil surface boundary temperature based on known 
ambient air temperature or in situ temperature measurements 
of the soil surface. 

2. Compute an equivalent temperature acting across a known 
thermal resistance to produce a net flux into the soil surface. 
With this method soil surface temperatures are computed rather 
than imposed. 

3. Mathematically model each of the individual thermal processes 
as precisely as possible and solve the resulting system of 
nonlinear equations simultaneously with equation (7). 

The method formulated in this report for simulating soil surface 

energy exchanges is a hybrid combination of alternatives 2 and 3. 


Individual thermal processes are modeled on a macroscopic scale, con- 
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verted to an equivalent temperature, and coupled to a variable turbulent 
heat transfer coefficient to simulate heat fluxes transmitted across the 
air-surface interface. This approach generates surface temperatures 


and eliminates the necessity of solving a system of nonlinear equations. 


3.2 Modeling the Individual Processes of Energy Exchange 
From a macroscopic viewpoint an energy balance at the air-ground 
interface can be approximated by the expression 


Qsoil * Qswr + Qniwr * Xconv ~ Qevap = 9 (14) 


where Q.,41 = heat flux into the soil column 


Qswr = Short wave radiation component 
Qniwr = net long wave radiation exchange term 
Qony = turbulent heat transfer component 


Qevap = evaporative heat loss term 
and where all heat fluxes are expressed in BTU/f£t2-day. 
Prior to discussing how equation (14) is incorporated into the variation- 
al analog for transient heat conduction, the methods employed to deter- 


mine each of the above energy components are briefly discussed. 


3.2.1 Short Wave Radiation 

Incident short wave radiation is determined by 

1. Applying CRREL nomographs (Gerdel et.al., 1954) for computing 

| available solar energy as a function of time and latitude. 

2. Modifying available solar energy with empirical relationships 
(Eagleson, 1970) which take into account altitude of cloud 


base and per cent cloud cover. 
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3. Selecting appropriate surface albedos from Wechsler and Glaser 


(1966). 


3.2.2 Long Wave Radiation Exchange 

As discussed in section 2.2.3, the net radiant heat emission term 
is nonlinear. If the surface emissivities of two radiating bodies are 
the same, a linearization technique discussed by Mettler (1966) and 
Pivovarov (1973) could be applied to provide 

Quiwe = 40°93, ($-..) (15) 
where all terms are as defined in section 2.2.3. 

Since the emissivity of a soil or snow surface is different than 
that of the atmosphere, it is easier to compute net long wave radiation 
exchange from the Stefan-Boltzmann law and apply net radiation heat 
transfer as a discrete flux across the boundary surface. An "apparent 
atmospheric emissivity" is assigned from empirical formulas (Eagleson, 
1970) which incorporate per cent cloud cover, cloud base altitude, and 
atmospheric water vapor pressure as influencing variables. 

Inaccuracies are introduced because previously computed boundary 
surface temperatures are used to determine net radiative heat transfer 
in the next time increment of the transient heat conduction problem. 
This is preferable to formulating an exact variational statement because 


computer solution by iteration is avoided. 


Swed Evaporative Heat Losses 


Evaporative heat transfer is not directly simulated but is estimated 
from data and applied as an outgoing flux across the surface boundary. 


The magnitude of evaporative flux for a given region is determined by 
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analyzing the potential evapotranspiration (PET) or actual evapotranspira- 
tion (AET) losses from data given by Patric and Black (1968). A time 
dependent evaporative heat flux function is constructed with a temporal 
shape closely approximating that of available solar energy. Its inte- 
grated area equals the selected PET or AET. In this report, tris 
desired to simulate the effect of evaporation from a saturated soil 
system; hence, PET is employed to compute an evaporative heat flux 
function. Evaporation during the winter months is assumed to be neglig- 


ible. 


3.2.4 Turbulent Heat Transfer 

Turbulent heat transfer is modeled according to the relationship 
- > 
where K = a turbulent heat transfer coefficient that is a 


function of wind velocity, surface roughness, and 
atmospheric stability 


Qonv * Kcureace air) (16) 


¢curface = surface temperature 
daiy = ambient air temperature 

The theoretical and experimental contributions of Scott (1964, with 

subsequent verification in 1969) are used to compute a turbulent heat 


transfer coefficient that is allowed to vary at the end of each time 


increment or when meteorological parameters change. 


3.2.5 Snow Cover Effects 

It is not desired to simulate the hydrologic morphology of sie! 
cover over the soil surface but to observe the gross insulating effect 
on a seasonal basis. In this respect heat conduction through an average 


winter snow cover is modeled using the previously discussed energy 






b AWN) Laat tne ae Ye mets £20b 60? 
 ferogen: Raith Gaktaamies bb amb ratwe gactt Hoi oy i reo 

a a ‘S21 -ygreno, talox othe bEeee Ro dads ant zesiicorie Mi 
7 Tig tt <axbiter ebm ol TRA oy TRE fose tox of tn if _ 
“ai Saxewwee: a wt Haksarereve ho ane wife om 
sartt peed oviteTégaye on etugibce oF bevciden shea i 
mod ot heewaas et eisndim temniw e472 yriteb rotten 
















5, itasnta ales ofa of gnibivsse bolotos ci totasinals sneit 1 
8H & “cen” ~-segaanal ‘+ ** 
eek Maid 2h0T SV es vetemet) Jad roles Be H. 
tue , Reenilgvot dddied .yttogioy. beim Yo noted 

Vitiidae uirtstyacess 
a 
OSUTATOGRDS OSHTSNA NW geet 
ea . ™ ae 
iL Sisteregves. 1k savidnas « nba 


- dtiw ,$001) J2092 to enaivatigines Bitnenticexe bm Das 















deed Ioeietbe 8 asughes o¢ bese oth [RU mt dotpgseaiteed theupel 
* enis doos Fo bea ads 3 fzey 09 Bewceils ax said tas lab Visss TO 2! 
Qycmio st totemene Ieotpouiexcetem wake to $e 







2708219 sevoD wane EAS a 
: Mone 2 To NG deere goliwond oly evofurts of Bovhaub son af 1h < , 
A p01 Nignivent seo%y ont svaeedo- os Tad soak” boa ads t4¥o seis | | 
rth as dauertiy nobtouhio jee st0ne5<e Shi  Lakeee ) 
“4 
‘a " “) i > gnianw Le tebon at -voves woo 
= 


29 


exchanges. Snow cover density is estimated from data presented by 
Biello (1969) and thermal conductivity from curves discussed by Mellor 
(1964). 

Minor sophistications are included to permit degradation of snow 
surface albedo when the ambient air feaperatiee rises above 32°F. 
Simulation of heat conduction through the snow cover is terminated 
when the computed snow surface temperature is greater than or equal to 


32°F. 


3.3 Finite Element Simulation of Surface Energy Exchange 

Having determined individual heat flux components, net energy 
exchange across a two-dimensional surface is simulated by coupling 
one-dimensional elements having a thermal conductivity equal to the 
turbulent heat transfer coefficient and having a volumetric heat 
capacity of zero to the two-dimensional elements representing the 


boundary surface. An equivalent temperature is calculated from 


: Qe Oke = Sevan 
047 tam * , (17) 


equivalent temperature 


= 
a 
1) 
4 
o 
©- 
0} 
ul 


= ambient air temperature 


a 


K = turbulent heat transfer coefficient 
and other terms remain as previously defined. The computed temperature 
from equation (17) is then imposed at the end of each one-dimensional 
element. t 
The imposition of equivalent temperatures at discrete nodes produces 


a net uniform flux (satisfying equation 14) across the two-dimensional 


boundary surface once they have been properly inserted into the system 
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matrix si Although substantiating mathematical details are 
omitted for the sake of brevity, note that if the equivalent tempera- 
ture is equal to the ambient air temperature, only turbulent heat 
transfer is simulated. Figure (6) is a schematic presentation of the 


concept. 


3.4 Simulation of Surface Energy Exchange under Fairbanks Field Conditions 
Numerical simulations of surface energy exchange were conducted 

for Fairbanks, Alaska summer and winter climatic conditions to permit 

evaluation of several solution techniques. Numerical instabilities 

encountered in modeling winter snow cover effects led to the development 

of the equivalent temperature approach for handling thermal processes 


at the air-surface interface. 


3.4.1 Surface Energy Simulation: 1-10 July 1972 

Using local meteorological data for Fairbanks, summer surface 
temperatures of exposed silt were simulated over a ten day period. Silt 
thermal properties were computed from the experimental data of Kerstens 
(1949). A clear day outgoing evaporative heat flux was determined from 
Fairbanks evapotranspiration loss data (Patric and Black, 1968) and 
further modified by a factor equal to the ratio of actual solar energy 
received to theoretical available clear day solar energy. Input data 
and results are shown in Table 1. Computed ground surface temperatures 


were approximately 8°F warmer than the ambient air temperature. 
t 


3.4.2 Surface Energy Simulation: 20-29 January 1973 
A ten day winter simulation of heat conduction in a soil column 


with an 18 inch snow cover was conducted to allow comparison of 
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NET FLUX INTO SOIL 
a oti sade JHA 





SYMBOLS. 
Pea = Equivalent Temperature 
ow-o = One-Dimensional Surface 
Energy Exchange Element 
K = Turbulent Heat Transfer Coefficient 
pe = Volumetric Heat Capacity=O 
Fig. 6 
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I. Meteorological Data . 
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Il. Simulation Results 
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Summer Simulation (1-10 July 1972) 
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measured and computed temperature profiles at the end of the 10 day 
period. Of secondary interest were the computed snow surface tempera- 
tures. Average Fairbanks January snow density and thermal properties 
were determined from values presented by Wheeler (1973). In this 
simulation, solar and long wave radiation fluxes were placed directly 
across the snow surface without employing an equivalent temperature 

and turbulent heat transfer coefficient. Solution instabilities caused 
by day to day fluctuations in meteorological parameters were encountered 
when all surface energy balance terms were simulated in such a manner. 
Available options to solve the stability problem were those of utiliz- 
ing very small simulation time increments or applying the equivalent 
temperature concept. The latter was chosen. 

Results of the winter simulation are shown in Table 2. Ten 
elements were used to simulate the snow cover and six elements to model 
a 45 inch soil column. A lower soil boundary temperature of 32°F was 
imposed since field observations indicated a mean January position of 
the groundwater table at a depth of 45 inches below the ground surface. 


Calculated and measured results were in good agreement. 
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IL. Simulation Results 
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It. Computed and Measured Soil 
Column Temperatures at the End 
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Chapter 4 


TWO-DIMENSIONAL PERMAFROST ISLAND THERMAL ANALYSIS 


4.1 Statement of the Problem 
A theoretical engineering analysis was conducted of the thermal 
feasibility of man-made permafrost islands. Two different problems 
were considered: 
1. The thermal degradation of an island drill site built on 
location during the winter months and completely frozen 
prior to spring breakup. 
2. The propagation of a winter freeze front through an island 
constructed from dredged seafloor sediments during the summer. 
Simulated island geometry was a right cylinder with a vertical 
height of 25 feet and radius of 24 feet immersed in 20 feet of water. 
Computer storage limitations prohibited investigation of larger island 
geometries. Two-dimensional finite element discretization of a vertical 
cross section of the island is shown in Figure (7). Symmetry in the 
transient heat conduction solution was assumed, thus, only the left 
half of the island is shown. 
Daily meteorological data from Barter Island, Alaska, was utilized 
to simulate surface boundary temperatures. Barter Island is located 
2 miles from the north shore of the Alaskan mainland in the Beaufort 
Sea. The climate is determined by the surrounding ocean environment 
which prevents extremely low winter temperatures and limits the warming 


effects of continuous daylight during the summer months. 
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4.2 First Thermal Analysis: Prefrozen Drill Site 

The first thermal analysis conducted assumed a prefrozen drill 
site built during the winter at an initial uniform temperature equal to 
the ambient air temperature (15°F). Thermal simulation was for the 
period 1 May 1970 through 31 August 1971. ~ 

The upper 4 feet of the island consisted of saturated gravel with 
a porosity of .45 and the remaining 21 feet were saturated silt with a 
porosity of .5. Soil thermal properties were determined from Kersten's 
(1949) report. Gravel in the upper part of the island that thawed 
during the summer was assumed to drain to a moisture content of 10 per 
cent of gravel unit dry weight. Draining would deter thawing, and 
gravel thermal properties were modified accordingly. Subsequent winter 
freezing was assumed to occur isothermally at 100 per cent soil satura- 


tion at 32°F. 


4.2.1 Applied Boundary Conditions 


To conserve use of computer time, upper surface boundary conditions 
were estimated by conducting a one-dimensional finite element surface 
energy exchange simulation. A vertical cross section of the island was 
employed. Barter Island meteorological data was provided in 24 hour 
increments, and daily surface temperatures were computed for use in the 
two-dimensional model. Outgoing summer evaporative heat flux was 
estimated from Figure (8). The presence of a 15 inch mean winter snow 
cover at a density of .33 g/cm? was simulated for the period 1 October 
1970 to 15 May 1971. Monthly averages of simulated energy terms are 


shown in Figure (9), where average net flux across the boundary surface 
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fe) MAY 30 JUNE 60° JULY 90 AUGUST !20 
TIME (DAYS) 


EVAPORATIVE HEAT FLUX (BTU/FT*® DAY) 


Fig. 8 


Estimated Barter Island, Alaska Evaporative Heat Flux 
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Fig. 9 
- Simulated Surface Energy Terms 
15" Winter Snow Cover Assumed 
Borter Island, Alaska 


(1 May 1970 -31 Aug. 1971) 
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is derived as a residual term. In Figure (10), simulated average 

* monthly surface temperatures and snow-ground interface temperatures 
are shown along with average measured ambient air temperatures for 
comparison purposes. 

The effect of a .02°F per foot geothermal gradient acting across 

a sea bottom of frozen silty clay was simulated and found to be neglig- 
ible. For this reason, the bottom of the two-dimensional island was 
considered to be nonconducting. 


A mean annual seawater temperature of 28.5°F was applied along 


the completely immersed vertical sides of the island. 


4.2.2 Simulation Results 

Two-dimensional computer simulation results are summarized in 
Figures (11) through (14). These figures show computed isotherms for 
the 1970-71 winter and following summer for the cases of mean seasonal 
and negligible winter snow cover. The second case assumes that snow is 
removed to promote maximum soil cooling. In both simulations, approxi- 
mately 2.5 feet of the surface gravel thawed by the end of the summers 
of 1970 and 1971. To determine porosity effects on thaw, a one- 
dimensional simulation predicted a 5 foot depth of thaw if gravel 


porosity was .30 instead of .50. 


4.3 seed Thermal Analysis: Dredged, Unfrozen Drill Site 

A second two-dimensional thermal analysis was undertaken to 
determine winter freezing front propagation in a dredged artificial 
island with the same geometry and soil structure as the previously 


discussed case. Natural extraction of heat was enhanced by assuming 
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Fig. 10 


Average Monthly Simulated Surface Temperatures 
and Measured Ambient Air Temperatures 
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Fig. I 
Computed Isotherms (°F) on | May 1971 
Simulation Time Span: 1! May 1970 - 31 Aug 1971 


15""Mean Annuol Winter Snow Cover Included. 
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Fig. 12 
Computed Isotherms (°F) on 31 May 1971 
Simulation Time Span: |! May 1970 - 31 Aug I97I 
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Fig. 13 
Computed Isotherms °F) on |! May !I97I 
Simulation Time Span: 1! May 1970 - 31 Aug 197! 


Winter Snow Cover Assumed to be Negligible. 
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Fig. 14 
Computed Isotherms (°F) on 31 Aug. 1971 
Simulation Time Span: 1! May 1970 - 31 Aug 1971 
Winter Snow Cover Assumed to be Negligible. 
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that no significant amount of snow would be allowed to accumulate on 
‘the surface of the island. 

Salinity of the seawater in the saturated silt and gravel-was 
assumed to be 15 parts per thousand which depressed the freezing point 
to 30.4°F. Oceanographic measurements (Hufford, 1974) indicate that 
salinities may very along the arctic coast from 3 to 20 parts per 


thousand in the upper 50 feet of water. 


4.3.1 Imposed Boundary Conditions 

As before, a one-dimensional thermal analysis of surface energy 
exchange was accomplished utilizing a saturated silt soil column to 
determine if enough silt would freeze over one winter to justify a two- 
dimensional thermal simulation. Results are shown in Figure (15) where 
depth of freezing front penetration is plotted against degree days of 
frost. It was judged that enough silt did freeze to consider a two- 
dimensional analysis worthwhile. 

It was also recognized that the imposition of a mean annual sea 
temperature of 28.5°F would not be completely representative of imposed 
winter boundary conditions along the vertical sides of the island. For 
this reason, the growth of sea ice was simulated and the temperature 
distribution in the ice was imposed along the upper submerged portion 
of the island considered to be in contact with the ice. Geothermal heat 


flux effects were again assumed to be negligible. 


4.3.2 Simulation Results 
Two-dimensional simulation results are shown in Figure (16). 


Approximately the upper 13 feet of the island froze during the first 
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Note: Computation of Cumulative 
Degree Days Based on 
a Datum of 32°F 
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Fig. 15 
- One-Dimensional Simulation of the Freezing 
of Soturoted Silt Using 
Barter Island Cimatological Dota oe 


(15 SEP. 1970 - 15 May 1971) 
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Fig. 16 
Two - Dimensional Simulation of the Freezing 
of oa Saturated Gravel and Silt Island 
Computed Isotherms (°F) on 15 May 1971 
Simulation Time Span: 1! Nov. 1970 - 15 May 1971 
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winter. To assist in predicting the elapsed time for the entire 
‘sland to freeze, a one-dimensional simulation incorporating geothermal 
heat flux was run with the results shown in Figure (17). It was 
estimated that the island would be completely frozen after 5.5 years. 
Extrapolation of results indicated an equilibrium frost penetration 
depth of 27 feet would be reached in a saturated silt soil structure 
with a 4 foot saturated gravel over-burden. 

A second one-dimensional simulation using saturated gravel pre- 
dicted that a gravel island of the same dimensions would be completely 
frozen in 3.5 years. Computed equilibrium freeze front penetration in 
gravel subjected to Barter Island climatological conditions was 48 feet 
after an elapsed time of 18 years. Simulation results are shown in 


Figure (18). 


4.4 Discussion of Results 

Both thermal simulations demonstrated that surface degradation of 
island drill sites due to summer thawing would be less than 3 feet, 
which is quite minimal. The primary explanation for this is that 
average Barter Island summer cloud cover is approximately eight-tenths, 
and consequently the magnitude of transmitted short wave radiation is 
low. As an example, average simulated transmitted solar heat flux, 
from Figure (9), during the month of June was only about 50 per cent 
of the available theoretical clear day solar energy. 

The convergence of computed isotherms (Figures (11) and (13)) near 


the air-water interface indicated that a large quantity of heat was 
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extracted from the surrounding water, transmitted through the island, 
and exhausted to the atmosphere .during the winter. 

Of considerable concern is the intrusion of heat along the vertical 
face of the island during the summer months as shown by the position of 
the 29°F isotherm at the end of the 1971 summer (Figure (14)). The 
danger of thawing can be eliminated by constructing prefrozen island 
drill sites with low salinity water. Dredged islands would not present 
a major problem since brine cells should migrate towards the freezing 
front during the winter and be rejected into the unfrozen silt. The 


high salinity silt would then become diluted during the summer months. 
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Chapter 5 


CONCLUSIONS AND RECOMMENDATIONS 


5.1 Conclusions 

The objective of this report was to formulate and develop a viable 
Tumerical model designed to solve arctic engineering problems concerning 
transient heat conduction accompanied by isothermal phase change. The 
first phase of the work involved comparison of the numerical model with 
several analytical solutions to phase transformation problems. Numeri- 
cal results and analytical solutions were found to be in good agreement. 

A second aspect of this research dealt with the simulation of 
surface boundary temperatures by employing meteorological parameters 
that are readily available for most geographic locations of interest. 
It must be emphasized that the theoretical and experimental results of 
numerous investigators were relied upon to evaluate the individual 
hydrologic energy exchange terms discussed in Chapter 3. Simulated 
energy terms and surface boundary temperatures were the results of 
simplified approximations to complex thermal processes. The numerical 
simulations of surface energy exchange under Fairbanks summer and winter 
climatic conditions assisted in assessing the simplifying assumptions 
employed. Of particular interest were the small discrepancies between 
simulated and measured soil column temperatures at the end of the 10 
day Fairbanks winter simulation, considering that ambient air tempera- 
ture fluctuated from -39°F to 16°F. 

The final phase of the work applied the thermal model to an 


engineering analysis of silt and gravel offshore drill sites. Based 
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on the thermal simulation results, the following conclusions are made 
subject to Barter Island, Alaska, climatological conditions: 
1. Under natural freezing conditions with negligible snow cover, 
equilibrium frost penetration depth for saturated silt with 
a 4 foot gravel overburden is approximately 27 feet and for 
Saturated gravel, 48 feet. 

2.. Summer thaw of saturated frozen gravel at a porosity of .45 
can be expected to be 3 feet and at a porosity of .30 can 
be expected to be 5 feet. 

3. Dredged island drill sites are thermally feasible in shallow 

water depths. 

Expanding on the above conclusions, the depth of water a dredged 
island could be built in would depend on a number of variables such as 
the magnitude of lateral sea ice forces and the manner in which the 
island is dredged. A possible method of constructing a dredged island 
is to place the silt or gravel in an expendable metal torus perforated 
to relieve excess pore water pressures as the island freezes. Use of 
a torus would prevent loss of material due to wave erosion and hasten 
freezing in contrast to dredged soils deposited at a small natural 
angle of repose. Since both silt and gravel freeze faster than seawater, 
the metal torus would not be destroyed unles the compressive strength 
of the frozen materials was exceeded or the frictional resistance of 
the entire structure was less than the maximum expected lateral force 
due to sea ice movement. 

Lateral resistance to sea ice could be enhanced by driving piles 


at the center of the torus and by using sea bed foundation design methods 
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_ currently employed on gravity structures in the North Sea. The 

t addition of an inclined rubble breakwater around the island would 
permit rafting of sea ice and reduce the horizontal component of ice 
forces. Furthermore, where drill sites are built in extremely shallow 
water, it becomes possible for seafloor sediments beneath the sites to 


freeze and increase resistance to lateral movement. 


5.2 Recommendations for Future Research 

Thermal analysis of permafrost island drill sites is no longer a 
problem of academic interest alone. Several drill sites have been 
constructed and occupied by Imperial Oil Limited of Canada in very 
shallow portions of the Mackenzie Delta (Riley, 1974), and plans are 
being developed to construct a dredged island drill site in 25 feet of 
water during the summer of 1974. As arctic offshore exploration 
progresses into deeper water, much larger ice forces will be encountered 
relative to the nearshore areas where ice mobility is somewhat limited. 
It will soon become necessary to develop methods of constructing in a 
relatively short time span massive structures capable of withstanding 
substantial ice sheet pressures. An economically feasible structure 
worthy of design consideration would be a grounded ice or ice-aggregate 
island. 

Before any large scale prototype ice island is built, a considerable 
amount of engineering research effort will be directed toward solving 
the thermal problem of ice eatetiorn: An optimum water thickness will 
be sought which will be easy to apply, will minimize thermal degradation 


of ice already produced, and when repeatedly applied, will yield the 
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desired thickness of ice within the available time constraints. Of 
secondary interest will be the problem of determining the volume of 
seafloor sediments capable of freezing beneath the base of an ice 
island. This is of value in computing resistance to horizontal forces 
and estimating the size of the island. 

It is recommended that the finite element method be considered as 
a research tool for these problems. With minor modifications to the 
numerical model formulated in this report, useful design solutions to 
the ice accretion problem could be generated. No modifications are 
necessary to solve the seafloor sediment freezing problem. In addition 
to the thermal problems, it is also recommended that laboratory tests 
be conducted to evaluate freezing front boundary shear strengths of 


Beaufort Sea sediments at probable future offshore drill sites. 








-([dalisyves ‘oe btw 207 to ete 





a 


Roemaiey 942. guiniars st to getdowy ett of fbiv savas: 
: ile te Go. 0e8 i? agreed ’ Jona? Fo -oldaqes ‘snqaitven 3 
™ 1033 detassicet 3 , tintzwahos at outey “to 28 ent PY 
m7 {2 eo asia aga gntanp | 


vhrgtunO3e% 


















MAS 22 eaoltasitibop tools tH .cowldesy eaedt toy get 


))@r envltuloe miest au. ate > mt bese tanred: istos 4 

‘eb enoltasitibosm of =. be2ex402; Bix . onan 
y 5. f , We 

inthis melacr fessor? 3 ; yf ines * avier oO 





23032 “(t0%es Gol san 1° : 3s aang? idan lsartot 
9 ke aitagiesta Seeiis yxubicned + + ovenulines. o3 
: ; ' i 


7 tegie 41h. atodei? ea q 1) otaemtsen owed 


1l. 


57 


REFERENCES 


Allen, D.N. and R.T. Severn, "The Application of Relaxation Methods 
to the Solution of Non-elliptic Partial Differential Equations: II. 


The Solidification of Liquids," Quarterly Journal of Mechanics and 
Applied Mathematics, Vol. 5, No. 4, 1952, pp. 447-454. 


Berggren, W.P., "Prediction of Temperature Distribution in Frozen 


Soils," Transactions of the American Geophysical Union, Pt. 3, 1943, 
pp. 71-77. 


Biello, M.A., "Relationships Between Climate and Regional Variations 
in Snow Cover Density in North America," U.S. Army Cold Regions 
Research and Engineering Laboratory Research Report 267, Hanover, 
N.H., 1969, 20 pp. 


Biot, M.A., Variational Principles in Heat Transfer, Oxford Mathema- 
tical Monographs, Clarendon Press, Oxford, 1970, 185 pp. 


Carslaw, H.S., and J.C. Jaeger, Conduction of Heat in Solids, Second 
Edition, Clarendon Press, Oxford, 1959, Chapter 11. 


Eagleson, P.S., Dynamic Hydrology, McGraw Hill, New York, 1970, 
Chapter 3. 


Gerdel, R.W., M. Diamond, and K. Walsh, "Nomographs for Computation 
of Radiation Heat Supply," Snow, Ice, and Permafrost Research 
Establishment, Corps of Engineers, U.S. Army, Research Paper 8, 
February, 1954, 5 pp. 


Hufford, G.L., "Dissolved Oxygen and Nutrients along the North 
Alaskan Shelf," Proceedings of the Symposium on Beaufort Sea 
Coastal and Shelf Research, Arctic Institute of North America, San 
Francisco, California, January, 1974. 


Jumikis, A.R., Thermal Soil Mechanics, Rutgers University Press, 
New Brunswick, New Jersey, 1966, Appendix 16. 


Kersten, M.S., "Laboratory Research for the Determination of the 
Thermal Properties of Soils," Final Report, Engineering Experiment 
Station, University of Minnesota, 1949, 226 pp. 


Keune, R. and P. Hoekstra, "Calculating the Amount of Unfrozen 
Water in Frozen Ground from Moisture Characteristic Curves," U.S. 
Army Cold Regions Research and Engineering Laboratory Special 
Report 114, Hanover, N.H., 1967, 7 pp. 





~7agrsl, . ite ee 
bist igs 









aoFIssUqeo) sot fefed A bee . bed) <A hake . 
Pianeta Sey ms ang har :2o2 .wort " vitge® 780 AOL e ee Se 


(8 tequt peer ore htc tt 22.8 (waeantand io 2g30o me fel y 
5 ag 2 reer «% 


AIyW Ons atite sanniteM hae nogveBheviczeta" vcd 

662 24OFuse® no mui doce? Ot Bo ep itéssiuyt." 2h ® 

find ,2oiweeA drsot 20 eFiditealwi7oth slimes): 2 oie ath 
CYS) tomes. .ninrottiao- gs 3 

(22909 yrterevigl Sgoysos ,esbaaityoh itok Lacedl . iA) 
ol hase aoel veers) wot cope weet 


o nolzantewrsde ado 2c? (oneeeed vtodetote' 2.8 nesepeyt 
tabat magi pxivewnsned ,ttecsa Ls * eling 20 £9)3caghed Lae 
@7 Ok, i ,Stoedors Fo ve cenaviad .2OReeee 


osial to Foe i) Riki THINS LE sceeseehi .4 gee A vom okt ; 7 
rti LIT oo er oo ues M tt oi) ucresT ai sone . . 
in tee, Tore’7od fxs gnd win doxvsoseRd anole bio? vee : 


OEs yh. tevores ast 





12. 


13. 


14. 


15. 


16. 


iy. 


18. 


19. 


20. 


21. 


22 


23. 


58 


Kolesnikov, A.G., "On the Theory of Ice Accretion on Sea-Surface," 
in Voprosy morskikh gidrologicheskikh prognozov (Problems of marine 
hydrological forecasts), Vol. 5, No. 12, edited by V.K. Agenorova 
and I.A. Benashvili, Leningrad, Gidrometerorologicheskoe Izdatelvo, 
1946, pp. 109-147, (in Russian). 


Koopmans, R.W. and R.D. Miller, "Soil Freezing and Soil Water 
Characteristic Curves," Soil Science Society American Proceedings 
30(6), 1966, pp. 680-685. 


Lukianov, V.S. and M.D. Golovko, "Calculation of the Depth of Freeze 
in Soils (in Russian), Bull. 23, Union of Transp. Constr., Moscow, 
1957. (English Translation, National Technical Information Service, 
Springfield, Virginia, 1957). 


Mellor, M., "Properties of Snow,'' U.S. Army Cold Regions Research 
and Engineering Laboratory Monograph III-Al, Hanover, N.H., 1964, 
105 pp. 


Mettler, A.J., and N.S. Stehle, "Ice Engineering - Analysis of the 
Growth of Sea Ice,'' U.S. Naval Civil Engineering Laboratory Technical 
Report R 497, Port Hueneme, California, November, 1966, p. 15. 


Myers, G.E., Analytical Methods in Conduction Heat Transfer, McGraw 
Hill, New York, 1971, Chapter 9. 


Nakano, Y. and J. Brown, "Mathematical Modeling and Validation of 
the Thermal Regimes in Tundra Soils, Barrow, Alaska," Arctic and 


Alpine Research, 4 (1), 1972, pp. 19-38. 


Patric, J.H. and P.E. Black, "Potential Evapotranspiration and 
Climate in Alaska by Thornwaite's Classification," U.S.D.A. Forest 
Service Research Paper PNW-71, Juneau, Alaska, 1968, 28 pp. 


Pivovarov, A.A., Thermal Conditions in Freezing Lakes and Rivers, 
Halsted Press, New York, 1973, p. 19. 


Portnov, I.G., "Exact Solution of Freezing Problems with Arbitrary 
Variation on Fixed Boundary," Soviet Physics - Doklady, Vol. 7, 
No. 3, 1962, pp. 186-189 (English translation). 


Riley, J.G., "How Imperial Built First Arctic Island," Petroleum 


Engineer International, 46 (1), 1974, pp. 25-28. 


Scott, R.F., "Heat Transfer at the Air-Ground Interface with Special 
Reference to Airfield Pavements," Technical Report No. 63, Arctic 
Construction and Frost Effects Laboratory, Waltham, MA, January, 
1964, 131 pp. 












Maan nae ~ stab seuaah-a: , 
: odizen wo emotdar?) 0dsas ys wy Ss 
7  RNSEOmE RA Pi ¥ “a hy rib Shaw 
,@elarsbsl aodzentdixo LoroTut sacs bi 


s 
. - ey 
«2B SRo2 007 Ee i] 
a a nr 


fer? Yo. ftqq0 ett zo 4702 fatuntie” gece Ae 
“pieago ,.ttened . genial sokidl eS. , Pht, (oen 
Soivede nolriemcza! Igutpdoel ’ Leakeektsy Oras. 
+ Oke 
\ oh i 


\ ; hs - q 
: in . y ; ; de vr44 : 
PRE .2 MLL) Is ~E IT Ante BOM: as 

o 


- ¥ 
ae 
of =~ . 5 i? AG. Test rae col C ‘ elites: 
patadsot rexeds) gnigscoryed [evil seuat Mi 
! am (mari BITC LAS Se a 2203 






- o 7 a) < 
Weta .antenewT geal) wabsouban) af fbodroM Ja ‘ata Peis’ 
ee a et a ye pet No . ope gh “1 - h 
Lb psigan) ) ayaa OY wel 

f. 


¢ baw, gableho! Lagtzarenray. .reodd 0 bee i 
- J . . 
Tha i pA wOrTe ‘oti S Otis 2) + omg i 
— 7 2s 2 oe sis! Vs <F | 1/ if: ‘ as 
E isntigengtaddsyve leiiostot™ ,Aoste .2.9 bite i wht 


4 : 1,2 1 noe cof4tecant) E16) T vet phe EA ae 
: a0 | A eoaai. ¢i\-*Vs 147 doteeee® oo 


- 


> dy . a 4 “4 7 3 Ae 7, 
E ae fSens7 4 SRO TEbsOD Letiadt( A ee) 


ri a) ges Rar ie (s2eth he 21 


. 7 4 rh ee ; — } 
* + - 
5 4 a< i» 4 : 
5 ‘Th t 
; 3 ; Ph » ‘ 


; J “4 ee 
+i ~ oh i 
@ bebseqal wit’ ..2.6 elke ae 
vd i¢ ea 4:9 aptpe t yo tA ey 
rb to tll =n ee a Lig .: } 
* 7 


Pek SP 


i Ad 
j te icaeifer’t 4eell" oe 17008 es ‘i 


¥ 





, iy col Nis 
> ss ’ » Ste 2 : 
0 / = a4 A sas 9 an ae “J 
i ; - DAB MOLI stra emo 7 en Pasi. 
' 


i wet 





24 s 


20 


26. 


27. 


28. 


29. 


59 


Scott, R.F.,''Predicted Depth of Freeze or Thaw in Soils by 
Climatological Analysis of Cumulative Heat Flow," U.S. Army Cold 
Regions Research and Engineering Laboratory Report 195, Hanover, 
N.H., 1969, 46 pp. 


Sears, W., An Introduction to Thermodynamics, the Kinetic Theory of 
Gases, and Statistical Mechanics, Addison-Wesley Publishing Company, 
Reading, Massachusetts, 1964, 374 pp. 


Stefan, J., "On the Theory of Ice Construction, Especially on Ice 
Construction in the Polar Sea," in Akademie der Wissenschaften, 
Mathematisch - Naturwissenschaftliche Klasse, Sitzungs - Berichte, 
pt. 2a, Vol, 98, Vienna, 1889, pp. 965-983 (In German). 


Wechsler, A.E., and P.E. Glaser, "Surface Characteristics Effect ~ 
on Thermal Regime, Phase I," U.S. Army Cold Regions Research and 
Engineering Laboratory Special Report 88, Hanover, N.H., 1966, 


26 pp. 


Weinstock, R., Calculus of Variations, McGraw Hill, New York, 1952, 
Chapter 7. 


Wheeler, J.A., "Simulation of Heat Transfer from a Warm Pipeline, 
Buried in Permafrost," Paper 27b, 74th National Meeting, American 
Institute of Chemical Engineers, New Orleans, Louisiana, 11-15 March, 
1973, 31 pp. 


64 
ws 
« 
> 
par 
_ nat 
- 
7(LLe 
7 7 
bo 
- 
‘ 
? 
’ 
ni 
7 


$3 


ont % rf tie Oe ht Lulde, 


wackt ,peealk edstldaddoetousivn vn 
wae) af) 88+ Ene «Ty sO88i ernetY., 





















wee 


Ce ee 
eitod 2. wi aft to rawr 
2.0 W wht roa eh 
oO! .gucmss viosereie 

? 
ro 


it 22 beth 


4 





eyo 


wh > freaks Pr’, tsa alot ods 


igen 1 
HW . ravens ba sroqal tahage aa 
A is new + air. dseizes zo 
THE c ToTQRRTT Oo 39 Kol 











othe a He 


ee tgs 





4. Appendix at 


“MATHEMATICAL DETAILS | sea Mae eae 7. 








i A ethane 


; DEIPATHG JADTTAMETTAM 


Variational Solution to the Heat Conduction Equation 


It is stated in Chapter 2 that, upon minimization, the generalized 


functional 


x = FfK, Go)? + KGS? + aces - c)ohaxay + Ladds 


(A.1) 
+ Un(o? - 26.6)ds + ofLegd> - egdeo)ds 
where c is that portion of the boundary surface where values of $ are 
not imposed. 

satisfies the heat conduction equation and the associated ioduness 
conditions of radiative, convective, and specified heat flux. The 
transition from equation (A.1) to the heat conduction equation will be 
demonstrated. 

Assume at a given instant in time ae'ts invariant and let the 
first term in equation (A.1) be represented by any arbitrary functional 
of the form 

SLE (XY 595 ox by)dxdy (A.2) 


Taking an arbitrary variation of f and its derivatives, produces 


e of of 
6x = FEL GSO) + (3g) * Ges poe + £a5eds 
(A.3) 
+ fh($sd - 4..64)ds + rac - €atn8d)ds 
Because the 6 and 2 operators commute and a minimum of the first 
variation is desired, the previous equation becomes 
of of 32 ee. 
= LflagS * 59, 5e(8*) * 5g, gy (6H) Iéxdy 
(A.4) 


+ £1a8d + h($8> - 4.89) + (e978 - £0550) ]ds = 0 
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The last two terms in the first integral of equation (A.4) may be 


integrated by applying Green's theorem of integral calculus such that 


IGE 5260) + ea gy hod iaady = “S{COL EGE * g Gel léxdy 
a 


(A.5) 





* foo axss + 1yse 5m Jds 


where 1x and ly are the direction cosines in the x and y directions 


respectively 


Upon substitution equation (A.5) into equation (A.4) we have 


=f Solas - GL - By Gig) Janey 
2x (A.6) 


+ [6o[1xe= + 1y2£ + gq + hb - deo) + O(Endt - cxht)]ds = 0 
30, * 3d) B 


Since 6 may be any arbitrary variation in R and ¢ is not imposed onc, 


application of the fundamental lemma of variational calculus requires 


that in R 

of 0 of 0 of . 

36 ~ 3x°39,) ~ dy%9y) = ° (A.7) 
and along c 

1x35 4 SE +q+h(¢ - $.) + o(epo* - Eafe) = 0 (A.8) 


By setting f equal to the first term of equation (A.1), it can be 
seen that we have an equivalent solution to the heat conduction equation 


if 6x is indeed a minimun. 







oe ie ae ms 
wevissessb < ine = ony tet euclaep foltosxih oily stn 4h 
> ce ae (an tlov magewn 9” 







, 
> epell. oN ae n) soto oon (2-A) soionmpe wot zuakzedue. 
¥ K , > nts ® 7 
Ca Seal Ame Pre WHR G, Qa: hen 2 i 


oe | “= 


z 

j : : 
a Ry ° , i 
: 7 


| “ta othe x os oof * 
te wes = egy * (A - 9)@ +p +~p4yl 4 spertheny a 









¥v _ 


7 





_ 
4 
~ 
a. 

: Ae ; 

“4 borogm! toa wz J burs QR At poiteivny ciastidip tie od nae het 


"at ; : : ry a ) 
. Senlype? aviualac tenotteivey be smo! Larnenahen? oie To wotseok 


y 1€ 
a,A) 5 % a, r i T& 
(a, A) ven? *a3) {6 oii 2 
‘ aa Cy Ay MG 2tipe FJ Y : am j 
matte, bien Seed arti c " 





inst struct Yous 


~ tee ‘, 


okey yey colves: one 


wong Ais cS here ior 


ee. temperate 2 my" be, apeckd! edt oesi 


ty WAL 


M Vi PRE oe ait Lt ave tor ae ot 


aE aie avert ace ty “prov i ding spec. 








© & elineqgh 


MAROON S3TUSO) 





B-1 


User Instructions 
a This finite element program solves one and two-dimensional 
transient heat conduction problems involving isothermal phase trans- 
formation. Boundary temperatures may be specified or internally genera- 
ted at the air-surface interface by providing specified meteorological 
parameters. One hundred and fifty elements and 100 nodes are permitted 
with the restriction that the numerical difference between any two 
numbered nodes common to an element not exceed 19; i.e., a maximum 
system bandwidth of 20 is permitted. All nodes must be consecutively 
numbered beginning with one and all elements must remain within the 
first quadrant of the Cartesian coordinate system. If one-dimensional 
elements are used, they are required to be vertically arranged. 

Energy simulation across a two-dimensional surface is accomplished 
by the use of one-dimensional elements as discussed in Chapter 3. The 
one-dimensional surface nodes and elements must be consecutively 
numbered from left to right beginning with one. The remaining two- 
dimensional elements and their respective nodes may be numbered in any 
desirable sequence. Do not specify imposed boundary temperatures on any 
of the one-dimensional nodes used for surface energy simulation. 
Equivalent temperatures are computed and imposed within the program 
itself. During data input assign a zero volumetric heat capacity and 
any fictitious conductivity to the one-dimensional elements. New 
conductivities are computed after meteorological data has been read in. 

Thermal units of BTU/FT/DAY and °F must be employed if the surface 
energy simulation mode of the program is selected. Any system of 
consistent units may be used for imposed boundary value problems with 


the exception that time be expressed in days. 
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Data Input 
FIRST CONTROL CARD (613, 5F10.2) 


A. 


cc. 


cc. 
Cc. 


cc. 


cc. 


cc. 


cc. 


cc. 


cc. 


cc. 


cc. 


1-3 


13-15 


16-18 


19-28 


29-38 


39-48 


49-58 


59-68 


1 for one-dimensional finite element simulation; 
2 for two-dimensional finite element simulation; 


3 for coupled one and two-dimensional elements for 


surface energy simulation. 

Number of elements. 

Number of specified temperature boundary nodes. 
Number of nodes not held at a specified temperature 
less number of nodes used for surface energy simulation. 
Total number of nodes. 

Number of one-dimensional elements used for surface 
energy simulation. 

Latent heat of transformation per unit weight of 
substance that is subject to phase transformation. 
Temperature at which isothermal phase transformation 
occurs. 

Time step size in DAYS or decimal fractions thereof. 
Desired elapsed time between nodal temperature print- 
outs, expressed in days. 


Total simulation duration in days. 


SECOND CONTROL CARD (15) 


cc. 


1-5 


Number of elapsed time segments before new boundary 


or meteorological data is to be read in. 
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THIRD CONTROL CARD (1F4.1); (Delete if surface energy simulation 
not desired) 


cc. 1-4 Average snow depth present. If no snow cover.is being 
simulated, enter a zero. 

FOURTH CONTROL CARD (1F5.2); (Delete if surface energy simulation 

not desired) 

cc. 1-5 Enter 1 if it is desired to reform system matrices in 
every time increment due to large fluctuations in 


meteorological parameters; otherwise enter a zero. 


FIRST DATA CARD (8F7.2) 
CC. 1-7 List initial temperature of nodes, excluding specified 

temperature boundary nodes, in node number sequence. 
CC. 8-56 Repeat above procedure. Use additional cards as 


necessary. 


SECOND DATA CARD (2F10.2) 
cc. 1-10 X coordinate of node. 
cc. 11-20 Y coordinate of node. 


CARDS MUST BE IN NODE NUMBER SEQUENCE! 


THIRD DATA CARD (313, 8F5S.1) 


cc. 1-3 First node number of element. 
CC. 4-6 Second node number of element. 
cC. 7-9 Third node number of element. Enter a zero if the 


element is one-dimensional. 
Node numbers of two-dimensional elements must be read 


in a counter clockwise direction around the element. 
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Cc. 10-14 Existing thermal conductivity of element. 


cc. 15-19 Thawed thermal conductivity of element. 


CC. 20-24 Frozen thermal conductivity of element. 
CC. 25-29 Existing volumetric heat capacity of element. 
CC. 30-34 Frozen volumetric heat capacity of element. 


CC. 35-39 Thawed volumetric heat capacity of element. 

CC. 40-44 Per cent moisture content of element. (Per cent of 
element unit weight). 

CC. 45-49 Element unit weight. 


ALL CARDS MUST BE IN ELEMENT NUMBER SEQUENCE! 


FOURTH DATA CARD (7(13,F7.2)) 

CC. 1-3 Node number. 

cc. 4-10 Specified boundary temperature. 

cc. 11-70 Repeat above sequence. 

DATA MUST BE PROVIDED IN NODE SEQUENCE! 

FIFTH DATA CARD (16,2X,7F8.2); (Delete if surface energy simulation 

not desired) 

CC. 1-6 Date (day, month, year), i.e., 071271. 

CC. 9-16 Ambient air temperature (°F). 

CC. 17-24 Wind velocity in miles per hour. 

CC. 25-32 Cloud cover in tenths. 

CC. 33-40 Cloud base altitude in 1000 feet. Enter 34 if 
ceiling is unlimited. 


CC. 41-48 Relative humidity in per cent. 
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CC. 49-56 Clear day solar heat flux (BTU/ft?-day). 


CC. 57-64 Outgoing evaporative heat flux (BTU/f£t?-day). 


J. REPEAT H and I based on the number of times data is to be updated 


or boundary conditions changed. 


Computer Program Listing 


Complete computer program listing and example output are provided 
at the end of this appendix. Printed energy balance terms of evaporative 
heat flux, net long wave radiation exchange, and convective heat flux 
are considered as heat losses from the air-surface interface when 


positive. 
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INTEGER DATE 


7025 FORMAT (15) 


7050 FORMAT(10X»* DATA ERRGR- UNSPECIFIED NODES+ BOUNDARY NODES DO 
1 NOT EQUAL TOTAL NUMBER OF NODES* ) <= raise 


7075 FORMAT (1F4.1) 
7100 FORMAT (8F7.2) 


7150 FORMAT(2F10.2) 

7200 
1250 
7275 

7300 

7350 
7400 
7450 

_ 7509 
7510 

7520 

ave, 325 
7530 

7550 

7600 

~~") COMPLETED® ) 


eS 


FORMAT(31378F5e1) 
FORMAT(T7(I3,F7.2)) | 
FORMAT(1692X97F8e2) 
FORMAT (10X,38HDATA ERROR-SYSTEM BAND WIDTH TOO LARGE) 
FORMAT(10X_,22HNEGATIVE AREA ELEMENT,I5) 
FORMAT(10Xy4HNODEs2X91592XeL3HFROZEN ON DAYs2XeF 702) 
FORMAT (10X~e4HNODE»2Xy1592XyL3HTHAWED ON DAYs2XeF 702) 
FORMAT(10X,"*COMPUTED NODAL TEMPERATURES ON DAY* y2XeF7e2//) 
FORMAT(7(2X_"*NODE'»3X_" TEMP.*)) i ascent ahaa are it 
FORMAT( 2X» 7( #—--—----—------")) 

FORMAT (//92X57( *----— eae 
FORMAT (2X,7( *-------——----- "),//) 
FORMAT(7(2X%914s3X2F 602) ) 

FORMAT (LOX, "SIMULATION OVER SPECIFIED TIME SPAN HAS BEEN 











8000 FORMAT (1LOXs3Xe"SNVP* 93Xs3Xe"QSOL? 93Xe2Xy"EVELX® 9 3Xe 2X "QLRAD', 


“8050 FORMAT(//,10X,"CGMPUTED SATURATION VAPOR PRESSURE AND ENERGY 


13X23X»2*QCONV* ,3X) 


IBALANCE TERMS CN DAY*,F7.2) 


8100 FORMAT (10Xs5F10.2) 
~ 8150 FORMAT(1LOXs"SNOW COVER HAS MELTED SNOW COVER SIMULATION TERMINATE 


1D ON DAY* sF7.29% e'e/910Xe*COMPUTED NODAL TEMPERATURES ARE * ) 


8175 FORMAT(1OXe"SIMULATION CALENCAR DATE IS(YRMODY) %32X,1Se//) 


~ 8200 FORMAT (1OX,*SOLUTION HAS BECCME UNSTABLE ENERGY SIMULATION HAS 


1BEEN TERMINATED® ) 


8225 FORMAT (10X,*COMPUTED EQUIVALENT TEMPERATURE®,» 1F6.2513X, "AMBIENT | 


LAIR TEMPERATURE", 1F6.2) 


8240 FORMAT{//,10Xe"THE FOLLOWING NODES ARE PARTIALLY FROZEN *) 


8250 FORMAT(10X,'NODE',5Xy"PERCENT OF ASSIGNED VOLUME FROZEN',5Xe*ASSIG __ 


1NED VOLUME * ) 


8275 FORMAT(1OXy14,19Xs1F10.2931Xy1F 522) 
8300 FORMAT(1OX,"DATA ERROR NO BOUNDARY NODE SPECIFIED. AT LEAST ONE _ 


1BOUNDARY CGNOITICN MUST BE READ IN AS DATA) 


8325 FORMAT(//,19Xs*NO NUDES UNDERGOING PHASE TRANSITION AT THIS TIME', 


~- 8350 FORMATI(1F5.2) 





“COMMON /BLK3/ P( 100,20) 


1//) OPES PSPCMI RNS SUNBATI [igtesyHiAs). oy 4 

COMMON /BLKL/ R100) 

COMMON /BLK2/ S(100,20) BEALE oe Sty SSS go 
COMMON /BLK4/ ST(190,20) 
COMMON /BLK5/ NOD(150,3) 
COMMON /BLK6/ X{(109) ay 
CCMMON /BLKT/ YC 100) 
COMMON /BLK8/ TK(150) 

“COMMON /BLK9/ ROWC(150) 
COMMON /BLKLO/ W(100,20) 
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_- 4s COMMON /BLKLL/ Z{1C0) é Le ao deebots SMe aes 
COMMON /BLK12/ NBC(59) 

COMMON /BLK13/ BC{50) : 

DIMENSION RWCF (LOO) »RWCT( LOO) sNKN{L00) »TO( 100) »G( 100) 

DIMENSION TXT( 159) sTKF{150)-RCWCT(150)sROWCF(L50)¢CMNOD(L50)_ 

LFLUID(150) »FNOD( 150), PChTR{1L50) pUNWAT(150) 

DIMENSICN VGL(100) ,FZHET( LOO) »NFREZ( 100) »WATE( 100) 

DIMENSICN A(3,3) »PEL(3,3)eXLCC(3),YLOC(3),SEL{3e3) 

READ(1L,7000) KODE» NEMS» LBC oe NK 9 MM, NSURF pHEAT p TFAZ »DT »DAY, END 

READ(1,7025) NTSEG 


COUNT=CAY/DT 











COUNTI=0. 
TIME = 0. 
Mam) hee fibAY=.0.. ag ain yar e = = <i a 
ITSEG = 0. 
SNOW = 0. 
J.  MODN: =. Osnye. +) lg am Sigh ace ee ie Po een eiwias 7 erae. 


IF «{NSURF) 200,200+100 
100 READ(1,7075) SNOW 
a READ( 128350) OPT 
200 CONTINUE 
IF (MM-(NK4+L8C) 2EQ. 0) GO TO 500 
WRITE(3,7C50) = — : : 
GO TO 3900 
500 CONTINUE 
aaa ea ~~ ¢—--=" Re alias aah aaa 
C—----READ INITAL NODAL TEMPERATURES 
C----- 
REAO(1,7100) € TO (I) » T= 1,NK )- 
00 1000 I=1,NK 
VOL{(I) = 0. 
Ls. WATECH) Ow. 
RWCF(I) 
RWCT(1) 
— IF (TO(L) «GT. TEAZ) GO TO 800 
NFREZ(1) =1 
: GO TO 1600 
7. 5 9s800 NFREZCE) = 2 
1000 CCNTINUE 
C----—- 
~~ €—==-READ FINITE ELEMENT MESH COORDINATES AND ELEMENT THERMAL 
C———--PARAMETERS 
C----- 
—, 7 REAOCL,TL5SO) CXOT) YET), T=1s¥M) 
READ(1s7200) ( (NOD(I,J) pJ=le3) eTKELT) sTKT( LT) e TKFC I) »ROWC(I), 
LROWCF( 1) ,ROWCT(1),PCWTR(I) UNWAT( I), I=L9NEMS) 
So s-~—CUp pe pe £005 1s kWENS 
 FNOD(I) = 0. 
1005 CMNOD(T) = 0. 
~~~ 1006 IF (LBC) 4600,%000,1007 
1007 NSERT = NSURF + 1 
IF ( NSERT «GT. LBC) GO TO 1113 
OD READ(197250,END =3800) (NBC(S),BC (J) »J=NSERT LBC) 
1113 IF { NSURF) 1LOC8,1C003,1111 
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\ 
tf ed vetoes oeuwew 


eae AED DONT V2) T= pi eh 


1l1l2 NBC(T) = Neh erNGEh ein GeVeRT Ser meet. 
C----—— 
2. G===-- DETERMINE UNKNOWN NODE POINTS __ 


Cc—--- 


1008 IF (TIME) 390041009,1355 
1009 IT =1 
NKN(I)=1 
1010 DO 1050 J=1,L8C 
LL=NBC(J) ; 
IF ( LL-NKN(I) .NE. 09 GO TO 1050 — 
NKN(I) = NKN(I) + 1 
1050 CCNTINUE 
TO6O) TH+ “Snowe ly“ tevkp ee eneesee hs eae oe a etree a 
IF ( I .GT. AK) GO TO 1070 
NKN(T)=NKN(I-1) + 1 
Tae “GO TO 1010 
1070 CONTINUE 
Cc--——- 
~“C—-—COMPUTE THE NUMBER OF UNSPECIFIED NODES COMMCN TO EACH ELEMENT 
C----- 
IF ( NSURF .LE. 0 ) GO TO 1300 
7 DO’ 1225°1=1,NSURF.. SEC Tere Re Lene | i UAE COR LANES een 
1225 CMNOD{(1) = 1. 
NSERT=NSURF#L 
GOSTO 1325000". Semi eee ie eo ch uso & nA coe a(S ee 
1300 NSERT=1 
1325 DO 1350 J=NSERT,NEMS 
#i,--- DGELSSOCKEL, Srereeneseen co Scene 
L= 0 
001350 M=1,LBC : 
IF (NBC(M)—-NOO(JyK) «FQ. 0) GO TO 1350 Sack STRSTR Tae 
IF ( NOD(JsK) .EQ. 0 ) GO TO 1350 
a a eee 
iwc i; .) Dente sotambecinGOl tO 41550...) (ose eo da mal i Le 
3 CMNOD(J) = CMNOD(J) + 1. 
1350 CONTINUE 
13SSCONTINUE 7 aoe, PRE RE CUNWN ouhi Minas Cea ho GR De. a i ana 
IF (NSURF) 1800,1800,1360 
1360 READ(1,7275,END=3800) DATE, TAMB,VEL yCLOUDyALT,RHy idee oie a5 eae 
“ers TF (VEL .GTe 1-) GO TO 1380 Seer ys. 
VELZ= 91. 
Cc ———— 
~ C—-—-ASSIGN APPROPRIATE THERMAL AND ROUGHNESS PARAMETERS 
c——— £ 
1380 IF { SNOW) 1385,1385,1395 
1385 EMISS = .90 
- ALBDO = .15 
GO TO 1430 
1395 EMISS = .$8 
IF(TAMB — 3220) 140021405,1405 
1400 ALBDO = .85 + elie 
an SGOMTOMIGUL «Me Goa Oe Meee a en ean 
"1405 IF (TOAY) 1410,1410,1415 
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Ss oft4 OT GO : «2 


: nbs ’ p864 BT GO ty 08. J 
geey UF O94 & rise 

— — nei 0) do 100d aS 
i + WN 


HEL. O00) 4°00! 
NAH Ay FIR UD 435 29 AAT CS PASO at erty 
UBEL Wi Val tet & 
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17 Bene Sls now Ad MANY FTALHUOATS 
OLE (VUES POE Coe TORE GREE 
we. = GEPMe 20€f : 
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~~ 1803 IF (TIME) 1805,1805,2020 


~ 3950 DO 2000 N=NDUM»NEMS _ 


1410 TOAY = TI 
1415 AGE = ABS 


1420 CONTINUE 


ME eS aah eT at See 
(TIME -— TDAY) 
ALBODO = 


Ct-—--- 


C—--—-COMPUTE CONVECTIVE HEAT TRANSFER COEFFICIENT. 


é——-—. 


1430 IF ( ABSUTAVE “- TAMB) .LT. 1.) GO TO 1440 


1435 
1436 CNVCT= 49.6¥*VEL*.38 © ’ 


IF (TAVE - TAMB) 14354144091445 
IF (SNOW) 14364143691437_ 





GO TO 1800 


_ 1437 CNVCT = 25.4#VEL*.38_ 


GO TO 1800 


1440 IF ( SNOW) 1441,1441,1442 


1441 


CNVCT = 86.6*VEL#.30 

GO TO 1800 Bare ars es 
CNVCT = 60.*VEL*.30 

GO TO 1800 

IF © SNOW) 14469144691447 _ 





1446 CNVCT = 120.*VEL*.22 


~~ 3900 CONTINUE © 


1447 CNVCT = 9T7.#VEL¥.222 


1800 REDO = 0. 


GO TO 1800 


GO TO 1800 


IF({OPT .LT. L.) GO TO 1810 
REDO = le 


1810 CONTINUE 


~ NROW = MM 


IF (NSURF) 180321803,1801 


1801 SUM = O. 


~ 00 1802 I = 1leNSURF 


1802 SUM = SUM + TO (1) 


TAVE = SUM/FLOAT(NSURF) 


WRITE (378200) 
GO TO 3900 


1805 NOUM=1 


C----- 
‘C——-—--CCMPUTE SYSTEM BAND WIDTH 


C---—- 


RCOk=Lic? He SVETER MATRICES 
IF { KODE .NE~ 1) GO TO 1900 
IF(KODE .€Q. 1) NCOL=2 

GO TO 2010 PRARLMEUAP: RODE DEL TY. A 
IF(KODE .£Q. 2) GO TO 1950 

NOUM= NSURF + L eViy oe 
DO 2000 I=1+3 

DO 20C0 J=153 ie eae oe 
NN = NCO(N,I) — NOD(N, J) * 2 


2000 IF « NCOL - NN .LT. O ) NCOL=NN 








62*EXP(—.0215*AGE) # .23¥*EXP(-.49*AGE) 


~ IF (€ ABS(TAVE - TAMB) .LT. 100.) GO TO 13803 
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_. 2010 thee ee 
NCOL «LE. 20 ) GO To 2020 
wale (3,7300) 


__GO TO 3900 
C---- 
C—---ZERO SYSTEM ARRAYS 
t---—- 
~ 2020 DO 2040 I=1yMM TAD eet yan pinata pbememeananne 
R(I) = 0. 


__00 2040 J=1,NCOL 
STC, J) 2°06 
S(I,J) = 0. 
P({IeJ) = 0. 
2040 W(I,J) = OO. 

IF (NSURF) 204592045,2041 

C---——- 

C—---COMPUTE INCIDENT SOLAR FLUX AND LONG WAVE ENERGY EXCHANGE 

C¢----- 

__ 2041 QSOL=(1.-ALBDO)*(1.—(.82-.024#ALT)*CLOUD/10.)#SRFLX 
TABS = (TAMB-32.)*(5./9.) + 273616 
TS=( TAVE-32.)*(52/92) + 273016 
VALUE = 373-16/TABS 

~SVPL = —7.902S8*(VALUE —-1.)# 5.02808*ALOGLO( VALUE) 

1—( 1.3816E-7) =( 10. #*(11.344%( 1.—-(1./VALUE)) )-Led 

24(8.1328E-3)¥*(10.**(-3.495149%(VALUE -1.)) 1.) 

“3 * 3.005715 TT ee ee 
SWVP= 10.**(SVPL ) 
VP=SWVP*RH/LOO. 

“TRANS= 1.0 — -O24*ALT 
ae _ ATMEM = 274 # .0125#CLOUD + .0049 *vP 

QLRAD=(1.—TRANS*CLOUD/10.) #(4.389E-7) *(EMESS#( TS##4. )-ATMEM* 

~ LI TABS##4.)) Zig, 


Cc---—- 
C—---CCNVERT RADIATION FLUXES TO AN EQUIVALENT TEMPERATURE 
as oe Sisets ao. ae eet or premier oe parimnew cn yet nedeoounalps 
TEQ = (QSOL — EVFLX — QLRAD) /(.336*CNVCT) 

EQTEM = TAMB + TEQ 
~~ "~2042 DO 2043 I=l»yNSURF 
BC(I) = EQTEM 

2043 TK(I) = .336*CNVCT 
Ct-—--—-- 
C—---CONSTRUCT THE SYSTEM MATRICES 
Cc---- 
“2045 DO 2600 M=1l,NEMS 
CALL CENTER (MyX@ARy YSAR eNSURF, KODE, DELTY» AM) 
IF ( KCDE .EQ. 1) GO TO 2550 
aapere ~~~ JF (M LE. NSURF) GU TO 2550 
CALL XIETA(M,XBAR» YBAReXLOC,YLOC) 
CALL AREA (XLOC,YLOC,AM) 
a BS BEL CAM 4GTs 023,60, TO, 2550 
WRITE(3,7350) M 
GO TO 3900 
2550 CONTINUE 
2050 IF ( M eLEe NSURF ) GO TO 2575 
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_.._ TF ( TIME .GT. 
DO 220C K = 193 
J=0 
OO 220C L = 1,L8C 
NCOM = NOD(M,K) 
IF ( NBC(L) — NCOM .NE. 0) GO TO 2075 
GO TO 2200 ; 

2075 IF ( NCCM .EQ. 0 ) GO TO 2200 

J=Jttl 


O-) GO.TO 2575 


2100 IF ( KCDE .EQ. 1 ) GO TO 2110 
' BETA = AM 
~—~ 60 TO 2120 

2110 BETA = ABS (DELTY)#*AM 
2120 DO 2125 I = L»NK 
— ‘FF -¢€ NKNCIY — NCGM .EQs 0 F GO TO 2130 
- 2125 CONTINUE 
2130 FLUID(M) = BETA*PCWIR(M) UNWATO(M) #(201) 
Rea e VOL({I) = VOL(I) + SETA/CMNOD(M) agai, 
: WATE(I) = WATE(I) + FLUIDUM) /CMNOO{M) 
zis RWCT(1T) = RWCT{I) + CROWCT(M) *8ETA)/CMNOO(M) 
—_ RWCF OL) = RWCE(E) + (ROWCF(M)*BETA)/CMNOD(M) ~ 
If{ TOCI) .GE. TFAZ) GO TO 2200 
FNOD(M) = FNOD(M) + le 
nes 22 OO SCONTINUE. cle ee merc hh aera nee 
2575 CONTINUE 
CALL SMALS(M,SEL sAMeg A gNSURF,KCDE,DELTY) 
CALL SMALP (M,AMePEL,NSURF,KODE,DELTY) 
CALL ADDIN (SELsPELeM,KODE) 
2600 CONTINUE : 
~~~" TF (TIME) 3900,2610,2E50 ~ 
2610 CONTINUE 
DO 2625 I=1:NK 
Tor S2ZHETCH) = WATECIISHEAT” 3 
RWCE(I) = RWCFII)/VOL(I) 
2625 RWCT(1) = RWCTCI)/VOLOI) 
fom nea eeoO CONTINUE |". pum mane ood 


C¢-—--- 





C—---INSERT SPECIFIED BOUNDARY TEMPERATURES AND MODIFY THE SePy AND 


ce. C—---R MATRICES ACCORDINGLY 


Ct---- 
CALL BNOIN{NROW,NCOL»LBC) 
SR a aa th FE ie) lalla 7: eee eae i EEC seo ) <2 oe 
C—--—--BEGIN TIME DOMAIN SOLUTION (CRANK NICOLSON) 
Ct---—- 





~~~ "p00 2700 f=1,NROW ae ery 

DO 2700 J=1,NCOL 

WCE, J)= SC1,J) * (2./0T) *P(TsJ) 
~~ 2700 Pll ygJ) = (2-/0TI (PKI eJ)) — SUIsJS) 


= 


C——--PRETRIANGULARIZE THE W MATRIX 
7 eee Cs = WurEt Agel: 4, eRe pt a Saat 


CALL PRESCL(NROW,NCOL) 
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“ a i; 
(IDK WORM SIZENS JAD VERE 


is. : aeRO s4ai5 


~ 3081 GII) = TFAZ 


~ C—-—==MODIFY ELEMENT THERMAL PARAMETERS WHEN NODAL PHASE CHANGE OCCURS ~~ 


.. 2800 ee eee 


~ 3084 CONTINUE 


~ .FZHETOL) = FZHET(I) + DELTQ  ~ 


3400 NCON = NKN(II) 


~ 3410 IF ( J.LE. NSURF) GO TO 3415— 
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CALL CCMB(TOeNROWeNCOL) 
DO 2900 I=l»yNROW 

___ 2900 Git) = 2.*R(T) + 201) 
CALL FINSOL (GyNROW,NCOL) 
TIME = TIME + OT 
COUNTI= CCUNTI + 1 
ITSEG = IISEG’+ L fee So ghy ys Thee Barca bed = Up hehe eee ers 
I=1 
GO TO 3084 


3082 TO (I) = G(I) 
I=I+1l 





IF €( IT .GT. NROW) GO TO 3689 
____—*ODIFF = G(I) — TFAZ 


“IF € DIFF) 3085,3C81,348T 
C---- 
C—-——-SIMULATE FREEZING PROCESS __ 
C¢---- 

3085 IF ( NFREZ(I) .EQ- 1 ) GO TO 3082 
DELQ’= RWCT(1)*VUL(1)*O1FF 
“FZHET(I) = FZHET(T) + DELQ 
IF ( FZHET(I) .GT. 0.) GO TO 3081 
NFREZ(I) = 1 

——~ GOI) = TRAZ + FZHETCID 7 ( RWCFCE)#VOL(I) 9 
WRITE(3,7400) NKN(I)y TIME 
FZHET (1) = WATECI)#HEAT 


GO TO 34000 -0000— am 





ee ee ee ee 


3487 IF ( NFREZ(I) NE. L) GO TO 3C€82 

3500 DELTQ = RWCF(1)*VOL(I)*OIFF*¥(-1.0) refi 
IF (FZHET(I) .GTe Of. ) GO TO 3081 

NFREZ(I) = 2 

G(T) = TFAZ — FZHET(I) / ( RWCTC(I)*VOL(I) ) 
WRITE(357450) NKN(I),TIME 

FZHET(1I) = WATECI) #HEAT 


DO 3415 J=1,NEMS 
DO 3415 K=1,3 
—NCHOS= NOD(J,K) 

IF ( NCHOS .£Q. NOON) GO TO 3410 

GO TO 3415 PEO BREE Os 
IF ( NFREZ(I) .€Q. 2 ) GO TO 3412 
FNOD(J) = FNOD(J) + Le 


3412 FNOD(J) = FNOD(J) — le 


C---- 


Cc---—— 


en: Neem _ ee a we) 
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pete __ 3413 uty = ( TK a J) #FNOD sel + { CMNCO(JI—-FNOD(J) )*TKT(J) 2)/CMNOD(S) 
ROWC(J) = ( ROWCF(J)*FNOD(J) + ( CMNOD(J)-FNOD{J) )¥*ROWCT(U) D7 
LCHNODE J) ' 
- __ 3415 CCNTINUE 
REOO = 1.0 aaa kn aie ter a eae ae 
GO TO 3082 
a __ 3689 IF (NSURF) 3694,3654, 3690 
c—— Mery crmmanner<acus 
C—---COMPUTE AVERAGE SURFACE TEMPERATURE d 
aa { 


ee — 


3690 SUM "Gz ee, 
DO 3691 I=1,NSURF 
__3691 SUM = SUM + TO(T) 
TAVE = SUM/FLOAT(NSURF) 
QCONV = CNVCT*.336*(TAVE — TAMB) ; 
IF(SNOW) 3694, 3654,3692 i 
3692 IF(TAVE — 32.) 3694736933693 fa? SESS eae Oa oe 
3693 SNOW = 0. 
WRITE(3,48150) TIME 
Se tS WMRETECs 15259) a Ee ee cat, ae) ore I eo eee one ee en 
WRITE(3,7550) € NKNOI)» TOCI) pI=1yNK) 
WRITE(347520) 
Mids cn -. 2 GO TO 35900 
3694 IF (TIME .GE. END) GO TO 3860 
IF (COUNTI .LT. COUNT) GO TO 3695 : 
= GOUNTI =50. 9. °° o alee SGI aa a ay gag Saar ae oT 

WRITE ( 3,7500) TIME 
WRITE(327520) 

~ WRITE(3,7510) : 
WRITE( 357520) i 
“WRITE(3s7550) ( NKNUIT)eTOUT)s I=LeNK) 
WRITE(3,7520) TPs Mg dl 
IF (NSURF .LT. L) GO TO 3695 

WRITE(328050) TIME 

“WRETEC OL SOAR Gao the lia. flee rae 
WRITE(3,8225) EQTEMsTAMB 

WRITE(3,7520) 

— WRITE (3 98000) — 
WRITE(347520) 
WRITE(3,8100) SWVPsQSOL,EVFLXsQLRAD,QCONV 
“WRITE(327530) 

3695 IFLITSEG «LT. NTSEG) GO TO 3700 
ITSEG = 0 
——~~~3700 IF (ITSEG) 371023710,3720 
3710 GO TO 1006 
3720 IF (RECO .GE. 1.) GO TO 1380 
mit GH MMe sUOMeCmEEIrr foc st) ao raar. 
3800 KRITE (347600) 
WREITE(397500) TIME 
i 2. © ONSURFSET<: 1).-GO: TO 3801 
WRITE(3,8175) DATE 
3801 WRITE(3,7520) _ 
ae WRITE(3e7510) 
WRITE(3,7520) 


—-}- ewan coows 
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EMM Bek gFIOT, (10904 | £00 Se: 
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QG06 OF 03 (sh «Tbe ; 

— WALLY (Ghee, CISte 
a a S70 . (272 a58) 
GMATLASTOR (enGeyO@39F 
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Ee 1 ree ay : (OOGGrE DBP yi. 
(OL2T +239 ee 
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a te 


(CCST +t ATI AW em 
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a , ‘ OFFER OLTEVOITE 132713 3 ae z 
. so: OF C8 OFTED 


Oe, Of O39 4,5 38. 6295) 41 -O8TE 


(o067 »t3 ou? 1 aw 

j gmM}7 COHET.ERST Day 
" (ee GT GO th} «Tj. WueM) 31 
STAG (ES 16.2927) Se 

40G2T. 29011 Pw 

wi a. Pay hy (Gset ef L9T 1 Aw 
(OER T LI 97 i kw 


min) oy ee yo mem + ach gene neg 


; i 
ok. ae ee ji up a “the — es ee eon 
, - . : . x a 
_ 7 t _ a 





~ TF (TQ(I) — TFAZ) 3810,3804,3810 





WRITE(3,7550) (€ NKNCI)eTOCI)s IT=lsNK) Pi A at ae, ae A 
WRITE(397520) 

MNTR= Q 

CO 3810 [=lsNK 


3804 MNTR = MNTR + 1 


__IF (MNTR NE. 1) GO TO 38C7_ 


WRITE(3,8240) 
WRITE(3,7525) 
WRITE (38250) 
WRITE(3,7520) j 


3807 CONTINUE 


CK= FZHETC(I)/(WATE(I) *HEAT) 

IF ( NFREZ(I) .EQ. 1) GO TO 3808 
PRCNT = (le — CK)¥*100. 

GO TO 3809 


3808 PRCNT = CK*1L00. 
3809 WRITE(3,8275) NKN(I)sPRCNTeVOL(I) 


3810 CONTINUE 


“Te (MNIR CETS-O] GO-10-3eLE Ose 
WRITE(3,-8325) 


3811 WRITE(3+7530) 


~ 3900 CONTINUE 


GO TO 5000 


_ 4000 WRITE(3,8300) 
5000 CONTINUE 


stop 
END 
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SUBROUTINE CENTER(M,XBAR »YBARyNSURF, KODE sDELTY, AM) 


Cc THIS SUBROUTINE COMPUTES THE CENTER COORDINATE OF TRIANGULAR. 
cl c ELEMENT M. 
__._. COMMON /BLK5/_NOD(15093) | 8g meee 
{ COMMON /BLK6/ -X(100) 
COMMON /BLK7/ Y(100) 


I = NOD(My1) _ 
J = NOD(M,2) 
IF ( KODE.EQ. 1) GO To 2 
IF (KODE .€Q. 2) GO TO1 
IF ( M .LE. NSURF) GO TO 2 


1 K = NOD(M,3) 


____ XBAR_= ( X(1) + X(J) + X(K) ) / 36 
~ YBAR = ( Y(I) + Y(J) + YCK) ) /3. 
‘ i GO T0 9 
2 DELTY = Y(J) - Y(T) 
AM = le 
IF ( KODE .NE. 3 ) GO TO 9- = 
AM = ( X(NSURF) — X(1)_) / FLOAT(NSURF) -__ a 
9 CONTINUE : 
i, RETURN 


——END 
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' COMMON /BLK7/ Y(100) 


SUBROUTINE XIETA (M,XBAR,YBAR »XLOCyYLOC) 

__ THIS SUBROUTINE COMPUTES THE LOCAL COORDINATES OF 
~ NODES ON TRIANGULAR ELEMENTS Me 

COMMON /BLK5/ NOD(150, 3) 5x: 

COMMON /BLK6/ X(100)_ er 
DIMENSION XLOC(3),YLOC(3) 
DO 100 L=1l+3__ 

I= NOD(MyL) 

XLOC(L)= X(I)-XBAR 


100 YLOC(L)= Y(I)— YBAR_ 


RETURN 2 
END 
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SUBROUTINE AREA (XLOC»YLOC,AM) 


THIS SUBROUTINE CALCULATES THE AREA OF TRIANGULAR 
__ ELEMENT M. een Mae) 
DIMENSION XLOC(3)+YLOC( 3) 
AM= (XLOC(2)*YLOC(3)—-XLOC(3)*YLOC(2)-—xLOC(1)*YLOC(3) 
__1 +XLOC(3)*YLOC(1) +XLOC(1)*YLOC(2)-xXLOC(2)*YLOC(1) )/2- 
RETURN ; : 
END 
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~~~ COMMON /BLK6/ X(100) 


—__ _____ SUBROUTINE SMALS(M,SEL,AMyA,NSURFyKODE,DELTY) __ 
c THIS SUBROUTINE COMPUTES THE SMALL S MATRIX 
C FOR TRIANGULAR ELEMENT M. 
COMMON /BLK5/ NOD(150,3) —_ Sons er 
COMMON /BLK7/ Y¥(100) 
_COMMON /BLK8/ TK(150) | 
DIMENSION SEL(393)9A(353) 
I = NOD(My1) 
J = NOD(M,2) Le a salen ene BA Seas 
IF ( KODE .EQ. 1 ) GO TO 3 = 
IF ( KODE .E0. 2) GO TO01 
__IF_( M .LE. NSURF) GO TO 3 _ 
1 K = NOD(M,3) 
Alls1)=X(J3)#Y¥ (K) =X (K)*Y (5) 
— ACLs 2)=X(K)*Y (CT) =X CT) FY (CK) 
A(193)=X (1) *#Y¥( 5) -X (5) *YV (1) 
Al2,1)= Y(S)-YCK) 
A(22)= Y(K)— YT) 


A(2,3) = Y(1) = YJ) 
A391) = X(K) — X(J) 

a A( 392) = X(I) — X(K) 7 
A(3,3) = X(J) - X(1) 
DO 2 J = 153 : 
DO 2 K = 1,3 


2 SEL(J9K) = Al2,J)#AC2,K)¥(TK(M)/CAM*S.) D 
1 + A(3,J)¥A(3,K)#( TKI(M)/(AM¥4.) ) 


GO TO 5 
eS tSsEUMS. = € TKUMEFAMY ZDELTYS 9°. 9 boo 
SEL(1,1) = ELMS*1l. 
SEL(152) = ELMS*(-1.) 
>. SEL(291) = ELMS*(-1.) 
SEL(2,2) = ELMS¥*l. 
DO 4 1 = 193 s 
meereeSELIC Isa). =: Os. ou ain RE a hn eee en Ei 
5 CONTINUE 
RETURN ecdnbas 
ENDO“. ot See 
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C 
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_____ SUBROUTINE SMALP(M,AM,PEL,NSURF »,KODE,DELTY ) | 2 
THIS SURROUTINE CONSTRUCTS A SYMMETRICAL P ELEMENT 
SUBSYSTEM MATRIX FOR THE MTH TRIANGULAR ELEMENT 

_COMMON /BLK9/ ROWC(150) 


DIMENSION PEL(3,3) 

IF ( KODE .EQ. 1 ) GO TO 2 
TE t RODE “60. '2!) Go. 70) 1.) 
IF ( M «LE. NSURF) GO TO 2 
EL = ( ROWC(M)*AM) / 12. 
PEL(1,1) = 2.*EL 
PEL (192)=1.*EL 
PEL (1,3)=1.*EL 
PEL(201)=1.3EL 2 
PEL (2,2)=2.¥*EL 
PEL(2,3)=1.*EL 
PEL (3,1)=1.%EL __ 
PEL (3,2)=1.¥*EL 
PEL (3,3)=2 EL 








GO SEO Ses a 2 ee 
2 ELP = ( AM*ROWC(M)*DELTY)/6~ 
PEL(1,1) = ELP#2. 
PEL(1+2) = ELP#*1. 
PEL(2+1) = ELP*1l. 
PEL (2.2) = ELP#2. 
_D0 3 T=1y,3_ 
3 PEL(I,3) = 0. 
5 CONTINUE 
RETURN 


END 
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SUBROUTINE ADDIN (SEL,PEL,M,KONE) 
COMMON /BLK2/ S(100,20) 
COMMON /BLK5/ NOD(150, 3) 

_____ COMMON /BLK3/ P(100,20) _ 

A DIMENSION SEL(3,3) 

DIMENSION PEL (3,3) 

SIPO (KODE EO!) 13) GO THY 200" 
IF ( KODE .EQ. 2 ) GO TO 1 

IF ( M .LE. NSURF) GO TO 2 

Se 2k DMNLOG) 32 1,3 

NR = NOD(MgJ) 

DO 100 K = 193 

NC = NOD(MsK) - NOD(M,J) + 1 

IF ( NC .LE. 0) GO TO 100 

S(NR»NC) = S(NR,NC) + SEL (JK) 

P(NRyNC) = P(NRyNC) + PEL(JeK) 
mElOOLCONTINUE |... a Sc 

GO To 5 

2 DO 200 J=152 
“NR = NOD(MyJ) 

DO 200 K = 192 
NC = NOD(M,K) — NOD(MyJ) + 1 EN 
~ TF ( NC .LE. 0) GO TO 200 
S( NR»NC) = S(NR»NC) + SEL( Je K) 
P(NR»NC) = P(NR,NC) + PEL (J9K) 
Te AZOORCONRINGE |... 1 aukat wap wee 
5 CONTINUE ; 
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COMMON /BLK1/ R(100 


_SUBROUTINE BNDIN (NROWsNCOL,LBC) _ 


) 


COMMON /BLK2/ S(100,20) 


COMMON /BLK3/ P(100 
COMMON /BLK12/ NBC( 


720) __ 
50) 


COMMON /BLK13/ BC(50) 


C- SET BOUNDARY CONDITIONS IN MATRI 


C- 


_00 3230 I=leLBC___ 
NN=NBC(T) 
NC=NCOL 
____3200_NR= NN=-NC +1 


TF(NR .LT. 1) GO TO 3210 
R(NR)= R(NR) — S(NR»NC)*BC(T) 


3210 NC=NC-1 
NC=1 
NR=NN 
3220 NR=NR41 


~ IF(NC «GT. 1) GO TO 3200 


IF {NR .GT. NROW) GO TO 3230 


NC=NC+1 


IF (NC «GT. NCOL) GO TO 3230 
R(NR)= R(NR)-S(NN»NC)*BC(I) 


GO TO 3220 


3230 CONTINUE — 
C- 
___C=REFORM MATRICES AND ELIMINATE EQUATIONS 
~C=AT BOUNDARY CONDITION NODES 


Cc 





___ 3280 _S(NR,NCOL)=0. 





DO 3300 [=1,LBC_ 
NN=NBC (I) —-I+1 
NROW=NROW-1 


__IF (NN .GE. NROW) GO TO 3250 


~ DO 3240 NR=NNyNROW 
R(NR)=R(NR +1) 
__DO 3240 NC=1,NCOL 
~S(NR »NC)=S (NR+1,NC) 
3240 P(NR,»NC)=P(NR+1yNC) 
___3250_NR=NN_ 
“N=2 
3260 NR=NR-1 
_IF(NR.LT. 1) GO TO 
N=N+1 


3300 


IF( N .GT. NCOL) GO TO 3280 


DO 3270 J=NyNCOL 
~ P(NRgJ—-1)=P(NRoJ) 
3270 S(NRyJ—-1)=S(NRo J) 
P(NRyNCOL )=0. _ 


S$ (NR yNCOL) = 0. 
GO TO 3260 


P(NRyNCOL )=0. 
3300 CONTINUE 
RETURN. 
END ti 
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_ SUBROUTINE PRESOL (NROW,NCOL) _ 


,_. C7THIS SUBROUTINE PRE-TRIANGULARIZES A SYMMETRIC 
i C-MATRIC IN BAND FORM FOR SOLUTION BY THE GAUSSIAN 
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C-ELIMINATION METHOD. FINAL SOLUTION IS BY SUBROUTINE FINSOL. 


COMMON /BLK4/ ST(100,20) 
COMMON /BLK1LO/ W(100,20) 
N=0 i 
100 N=N+1 
_STUNyNCOL) = W(Ny1) _ 


a 
C-REDOUCE PIVOT EQUATION 
c- ba 3: eles ie Rig ote ror 
IF(N-NROW) 150,300,150 
150 DO 200 K=2,NCOL 
___STOUNsK-1)=W(N,K) 
200 WINsK)= WINsK)/W(Ng Ll) 


C- 
___C-REDUCE REMAINING EQUATIONS 
Cc- 
DO 260 L=2+NCOL 
od AOA re 296 
ioe IF(NROW-I) 260,240,240 
240 J=0 
DOSES k= NCOs 
J=J+1 


250 WIT sJ)=W(I eJI-ST(NeL—-1) *W(N eK) 
260 CONTINUE 





GO TO 100 ae 
300 RETURN 
cay END ‘ceiving er eeniernen eae 





—————  . 
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a (------------— ------------------------------------------------- - 
__C-_THIS SUBROUTINE MULTIPLIES SYMMETRIC MATRIX P__ we 
C-TO MATRIX Y AND STORES THE RESULT IN Ze 71 7 es i 
Cc wee ewe ow ww wn a a ow oe ww we ww oe ww we ow ow wr we wwe oe owe ow ww wwe owe oe ow we we ww ww ew owe 
nal ger : 





“COMMON /BLK3/. P(100,20) 
COMMON /BLK11/ Z(100) 
_DIMENSION Y(100) 
DO 200 I=lyNROW 
‘ Z(T)=VCT)*P CT, 1) 
DO 200 K_= 2,yNCOL_ 
L=I-K+1 
IF(L.LT.1) GO TO 100 
_ZU1)=Z(1) + YCL)*P (LK) , ; 
“100 N=I+K-1 * MiNi een ea 
IF(N.GT.NROW) GO TO 200 
Z(T)=Z(1)_ +Y¥(N)#P CT) K) 


| 
| 











2 200 CONTINUE 

2 RETURN 
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ae C-THIS SUBROUTINE SOLVES A SET OF LINEAR SIMULTANEOUS 
,-- CTEQUATIONS WHOSE COEFFICIENT MATRIX HAS BEEN PRE-TRIANGULARIZED 
}  C=BY THE GAUSSIAN ELIMINATIION METHOD. THE SYSTEM MATRIX 
t  C=-IS IN BAND FORM AND IS SYMMETRICAL. SOLUTION IS PLACED 
__C-IN THE LOAD MATRIX. ; 


COMMON /BLK4/ $T(100120) 
COMMON /BLK10/ W( 100,20) 
DIMENSION S$S(100) 


C- 
~~ €=REDUCE LOAD VARTABLES 
C- 


N=0 
———~ 1007 N=NFT 

SS(N) = SS(N)/ST(NyNCOL) 

IF (N-NROW) 110,200,110 





‘110 CONTINUE 
DO 130 K=2,NCOL 
L=N+K-1 
~~ TF(NROW-L) 130,120,120 
120 SS(L)=SS(L)—ST(N,K-1)#SS(N) 
130 CONTINUE 
—~~~60° TO 100 c 
C- 
__C-BACK SUBSTITUTION _ a a , 
C- 
200 N=NROW 
300 N=N-1 


'” TFIN) 350,500,350 
350 DO 400 K=2~*NCOL 
L=N+K-1 
IF(NROW-L) 400,370,370 
370 SS(N)= SS(N) — W(NeK)*SS(L) 





i 400 CONTINUE 
= seen 300 >. 5 
500 RETURN 
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